Automated Efficient Estimation using Monte Carlo Efficient Influence Functions

Authors: raj  agrawal
sam  witty
andy  zane
eli  bingham
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 09:07:48

Abstract: Many practical problems involve estimating low dimensional statistical quantities with high-dimensional models and datasets. Several approaches address these estimation tasks based on the theory of influence functions, such as debiased/double ML or targeted minimum loss estimation. We introduce \textit{Monte Carlo Efficient Influence Functions} (MC-EIF), a fully automated technique for approximating efficient influence functions that integrates seamlessly with existing differentiable probabilistic programming systems. MC-EIF automates efficient statistical estimation for a broad class of models and functionals that previously required rigorous custom analysis. We prove that MC-EIF is consistent, and that estimators using MC-EIF achieve optimal $\sqrt{N}$ convergence rates. We show empirically that estimators using MC-EIF are at parity with estimators using analytic EIFs. Finally, we present a novel capstone example using MC-EIF for optimal portfolio selection.

Main paper: 
 introduction: 
 Over the past several decades, there has been remarkable progress on robust and efficient statistical estimation, especially for high dimensional problems. A particularly compelling class of such methods are built on a foundation of efficient influence functions (EIF), i.e., functional derivatives in the space of probability distributions [Ken22] . These methods have been particularly fruitful in causal inference applications, where estimating quantities such as the average treatment effect require modeling high-dimensional nuisance parameters relating confounders to treatment and outcome variables. Intuitively, these methods focus finite statistical resources on quantities that matter, and not on nuisance parameters that only indirectly inform the statistical quantities we wish to estimate.Despite their successes, estimation methods based on the EIF have lagged behind the kinds of automation that machine learning practitioners have grown accustomed to, instead requiring complex manual derivation on a case-by-case basis. This is contrasted with the generality of automatic differentiation (AD) systems [BPRS18] and probabilistic programming languages (PPLs) such as Pyro [BCJ + 19] or Gen [CTSLM19] , which automate numerical computations for probabilistic inference. EIF-based estimators have historically eluded this level of automation and generality, because exact recovery of the EIF requires solving high-dimensional integral equations.Contributions. We introduce Monte Carlo Efficient Influence Functions (MC-EIF), a general and automated technique for numerically computing EIFs using only quantities that are already available from existing AD and PPL systems. Our key insight is that EIFs can be expressed equivalently as a product of (i) the gradient of the functional, (ii) the inverse Fisher information matrix, and (iii) the gradient of the log-likelihood, as shown in Theorem 3.4 in Section 3. In Section 4, we show how MC-EIF can be used to automatically construct a variety of efficient estimators for a broad class of models and functionals, avoiding the need for complex manual and error-prone derivations.In summary, we show that: (i) MC-EIF provides accurate estimates of the true EIF, enabling efficient estimation, and (ii) MC-EIF is very general, applying to many functionals and models that can be written as probabilistic programs. In Section 3, we introduce MC-EIF and provide a non-asymptotic error bound on the quality of our approximation. We show how estimators using MC-EIF achieve the same asymptotic guarantees as using analytic EIFs in Section 4. In Section 5, we show empirically that MC-EIF produces more accurate estimates of the EIF than existing automated approaches, and using MC-EIF as a drop-in replacement for the analytic EIF does not degrade estimation accuracy in a variety of benchmarks, including a novel capstone on optimal portfolio allocation.Related Work. Influence function-based estimators have a rich history in the statistics and machine learning literature [BKB + 93, Tsi06, Ken22, HDDOV22]. Despite their effectiveness, these methods have historically required custom and complex mathematical analysis for specific combinations of models and functional. Even an incomplete recent survey of the influence function-based estimation literature yields a large collection of complex scenario-specific research. For targeted minimum loss estimation (TMLE) [VDLR06] ; cluster-randomized trials [BvdLA + 23], continuous time-dependent interventions [RGvdL22] , mixed experimental and observational data [DTA + 22], mediation analysis with longitudinal data [WvdLP + 23], subgroup treatment effect estimation [WPvdL + 23], survival and competing risks analysis [RvdL24] , continuous time-to-event outcomes [REvdL23] , and variable importance measures for effect estimates [LHvdL23] . For double/debiased machine learning [CCD + 18]; difference-in-differences [Cha20] , instrumental variable designs [JTB21a] , and mediation analysis [FHL + 22].. Importantly, our work does not introduce novel efficient estimators; instead it aims to lower the mathematical burden for practitioners who wish to use existing influence function-based efficient estimator templates (see Section 4) with custom models and/or functionals.Our work is not the first to attempt to automate and generalize computations for efficient statistical estimation. Perhaps the closest technique we are aware of is approximating the influence function using finite differences on kernel-smoothed empirical distributions [FQWD15, CLvdL19, JWZ22a] . We provide a thorough comparison with this method in Section 5, demonstrating how MC-EIF automates and scales better to high dimensional problems, exactly the settings where efficient estimation is most useful. Recent work has made progress towards general efficient estimators, but still impose strong restrictions on models and/or functionals. DML-ID [JTB21b] extends double machine learning to nonparametric causal graphs and marginal density under intervention functionals. Similarly, the kernel debiased plugin estimator [CGMM23] implements a version of TMLE that bypasses the influence function computations for models defined in a RKHS. Finally, other recent work [CNS22, FS23, CNQMS21, CNQMS22] approximates the efficient influence function for generalized method-of-moment estimators.Finally, we note that there are a number of intriguing connections between three related but distinct mathematical objects: the efficient influence function in semiparametric statistics [Tsi06] , the natural gradient in information geometry [Ama16] , and the so-called empirical influence function in robust statistics and machine learning [Law86] . A comprehensive review of these connections is beyond the scope of this paper, and we focus here on two particularly important points for contextualizing our work. First, we emphasize that the efficient and empirical influence functions are not equivalent: the efficient influence function is a fundamental mathematical object in semiparametric statistical theory which quantifies the effect of perturbing a functional in an arbitrary direction in the space of probability measures, while the empirical influence function is a distinct and somewhat more specialized [BNL + 22] object quantifying the effect of perturbing individual training points in a parametric statistical model. More specifically, in the parametric setting the efficient influence function is defined as shown in Theorem 3.4 in terms of the Fisher information matrix [Tsi06] , whereas the empirical influence function is defined in terms of the Hessian of a model at the training points, two quantities with very different mathematical and statistical properties that are not straightforwardly interchangeable [KHB19] . Second, we note that existing algorithms for computing natural gradients and empirical influence functions cannot immediately be adapted to efficient estimation. Specifically, our algorithm described in Section 3 for Monte Carlo approximation of the efficient influence function is similar in structure to some previous algorithms developed for estimating the natural gradient [TR19] and empirical influence function [KL17, GSL + 19]. This would seem to suggest porting other methods that compute more heavily biased approximations of the natural gradi-ent [GLB + 18] and empirical influence function to computing the efficient influence function, as some of these methods are known to scale to even the largest neural network models deployed in practice today [GBA + 23]. However, these approximations are not directly applicable in our setting because provably efficient estimation is only known to be possible with tight control over any approximation error introduced in computing the efficient influence function, as discussed in Section 4 below and in [JWZ22a] . Relaxing these restrictions to enable similarly scalable variations on the basic MC-EIF framework of Section 3 is an important direction for future work. 
problem statement: 
 General Problem. We consider the estimation of some estimand θ * ∈ R L , where L ∈ N denotes the dimension of the target quantity. Typically, we can express θ * = Ψ(P * ) for some known functional Ψ, where Ψ maps a probability distribution to a vector in R L , and P * (x) denotes the true data-generating distribution over some vector of observables x ∈ R D , D ∈ N. Many estimation tasks involve highdimensional nuisance parameters, or quantities of no immediate value to the analyst. For example, to estimate the average treatment effect, one might need to adjust for high-dimensional confounders.Semiparametric Solution. Semiparametric statistics provides a mathematical framework for optimally estimating θ * in the presence of potentially complex, high-dimensional nuisance parameters. A standard way to estimate θ * is with the plug-in approach; construct an estimate P of P * and report θ = ψ( P). Unfortunately, the plug-in approach can lead to provably sub-optimal estimates of θ * due to poor estimates of P * [Tsi06, CCD + 18, FS23]. Instead, a general recipe for efficiently estimating θ * from finite data {x n } N n=1 , wherex n iid ∼ P * (x) for n = 1, • • • , N, is given by the following three-steps: (i) use N/2 samples to construct an initial estimate P of P * , (ii) compute the influence function (to be defined shortly) of Ψ at the estimate P, and (iii) evaluate the influence function at the held out N/2 datapoints to derive a corrected estimate. 1 In Section 4, we elaborate on how influence functions are used to construct several popular efficient estimators.Influence Functions. A central premise of this paper is that to automate efficient estimation, it suffices to automate the computation of influence functions, which can be thought of as gradients in function space. We make this precise below.Definition 2.1. (Gateaux derivative) Consider the ϵ-perturbed probability distribution P ϵ := (1ϵ)P + ϵQ = P + ϵ(Q -P), where Q is some probability distribution. Ψ is Gateaux differentiable at Q if the following limit exists:d dϵ Ψ(P ϵ ) ϵ=0 = lim ϵ→0 Ψ(P ϵ ) -Ψ(P) ϵ .The Gateaux derivative can be viewed as a generalization of the directional derivative from ordinary calculus; it characterizes how much a functional changes at a point P in the direction Q -P.Definition 2.2. (Influence function) Suppose there exists a square integrable function φ ∈ L 2 (P)such that d dϵ Ψ(P ϵ ) ϵ=0 = ⟨φ, q -p⟩ L 2 = x∈R D φ(x)(q(x) -p(x))dx for all Q ∈ M and E x∼P [φ(x)] = 0, where M denotes some space of probability distributions and p(•) and q(•) are the density functions for Q and P, respectively. Then, φ is called an influence function for Ψ at P.An influence function is a re-centered "functional gradient" in L 2 (P): just as the Euclidean inner product between the gradient of a function and a vector yields the directional derivative in ordinary differential calculus, the L 2 (P) inner product between the influence function and perturbation Q -P yields the Gateaux directional derivative. Influence functions, however, are not always unique [Tsi06, Ken16] -some may lead to higher asymptotic variance estimators than others. The optimal influence function minimizes asymptotic variance, and is called the efficient influence function (EIF). When the EIF exists, it is P almost everywhere unique, and found through a Hilbert space projection onto what is known as the nuisance tangent space. We defer details to [Tsi06] and [Ken16] .As the influence function in Definition 2.2 is defined implicitly as a solution to an infinite set of integral constraints over M, it is often hard to find. Entire papers have been written to analytically derive influence functions; see, for example, the papers listed in Section 1. For even experts in machine learning and statistics, such derivations are out-of-reach, time consuming, and error prone.3 Monte Carlo Efficient Influence Function A key question we address is whether using a high-dimensional parametric model simplifies the process of solving Definition 2.2. We show that it does below.Notation. We let ϕ ∈ Φ ⊂ R p denote a finite-dimensional parameter specifying a distribution on the observed random variables x ∈ R D for p < ∞, p ∈ N. P ϕ (x) corresponds to a distribution in this space, and P ϕ * (x) the true distribution, or the one closest to the true data-generating distribution in Kullback-Leibler distance. We let ψ(ϕ) denote a function R p → R L that equals the evaluation of the functional Ψ(P ϕ ) for all ϕ ∈ Φ. Under mild differentiability assumptions, we provide the analytic formula for the EIF in Theorem 3.4.The first assumption states that the density of P ϕ (x) is continuous and differentiable with respect to ϕ, a condition satisfied by many parametric model families. For example, the univariate Gaussian density 1 √ 2π exp -0.5(xϕ) 2 is a continuous and differentiable function of its mean, ϕ ∈ R. Assumption 3.1. ∀x ∈ R D , the map ϕ → P ϕ (x) is continuous and differentiable with respect to ϕ.The next assumption is also satisfied for many functionals. For example, consider the mean functionalΨ(P ϕ ) = E x∼P ϕ [x].Continuing with the univariate Gaussian example from above, where the mean is unknown, we have ψ(ϕ) = ϕ, which is a continuous and differentiable function of ϕ. Assumption 3.2. ψ(ϕ) is a continuous and differentiable function of ϕ.The last assumption requires that the Fisher information matrix be invertible, which is necessary for ϕ to be identifiable [Tsi06] . Assumption 3.3. Fisher information I(ϕ) := E x∼P ϕ (x) [∇ ϕ log P ϕ (x)∇ ϕ log P ϕ (x) T ] is invertible. Theorem 3.4. (Theorem 3.5 in [Tsi06] ) Suppose Assumption 3.1, Assumption 3.2, and Assumption 3.3 hold. Then, the efficient influence function φ ϕ (x) at ϕ evaluated at the point x ∈ R D equals[∇ ϕ ψ(ϕ)] T I(ϕ) -1 ∇ ϕ log P ϕ (x).(1)While Equation 1 has been around for many decades, it has mainly been used as a theoretical tool for mathematical statisticians. In particular, Equation 1 is typically evaluated at the true data generating parameter ϕ * to characterize the theoretical asymptotic variance of an estimator. In other instances, it is used to derive approximate confidence intervals; see, for example, Chapter 3 in [Tsi06] . In the following Sections, we discuss how Equation 1 provides a key ingredient in automating efficient estimation in high-dimensional parametric models. 
numerically approximating the eif: 
 Given a model P ϕ (•) and functional Ψ(•), we seek to automatically compute Equation 1. Our Monte Carlo efficient influence function (MC-EIF) estimator achieves this automation by replacing ψ(ϕ) and I(ϕ), which are typically unknown, with stochastic approximations ψM (ϕ) and ÎM (ϕ) computed from M ∈ N Monte Carlo samples:EQUATIONHere, we show that Equation 2 leads to an automated and accurate approach to numerically computing EIFs using only quantities provided by existing AD and PPL systems.Approximating ÎM (ϕ) -1 ∇ ϕ log P ϕ (x). We draw x m iid ∼ P ϕ (x), 1 ≤ m ≤ M for M ∈ N, and letEQUATIONA naive approach for computing ÎM (ϕ) -1 ∇ ϕ log P ϕ (x) is calculating the full p × p matrix in Equation 3, inverting it, and then taking its product with the score vector ∇ ϕ log P ϕ (x m ) T ∈ R p computed from AD. This naive approach takes O(M p2 + p 3 ) time and O(p 2 ) memory which might be too expensive for large p. Instead, we exploit AD and numerical linear algebra techniques to avoid explicitly storing and inverting the approximate Fisher information matrix, similar to [KL17] . Suppose that we have a black-box method to compute Fisher vector products ÎM (ϕ)v for arbitrary vectors v ∈ R p . Then, we could use the conjugate gradient algorithm to iteratively find ÎM (ϕ)∇ -1 ϕ log P ϕ (x), where the cost of each conjugate gradient step is determined by the cost to compute ÎM (ϕ)v. While the number of conjugate gradient steps needs to be p for an exact inverse, often far fewer iterations are required for a close approximate solution [WPG + 19]. To make computing ÎM (ϕ)v efficient, we collect the M simulated datapoints in the matrix X M ∈ R M ×D and letEQUATIONwhere J M = ∇ ϕ log P ϕ (X M ) ∈ R M ×p is the Jacobian matrix. We use Pearlmutter's trick to avoid computing the entire Jacobian matrix [Pea94] . In particular, this method allows us to compute the Jacobian vector product v M = [J M v] ∈ R M in time proportional to a single evaluation of log P ϕ (X) and O(M + p) memory. Similarly, we use the vector Jacobian product to compute J T M v M . Approximating ∇ ϕ ψM (ϕ). Robust estimation with MC-EIF does not require exact gradients. Instead, it only requires a sequence of gradient estimators {∇ ϕ ψm (ϕ)} ∞ m=1 of ∇ ϕ ψ(ϕ) whose error can be bounded above by some ∆ m > 0, where the M th iterate ∇ ϕ ψM (ϕ) is used in Equation 2. 2 Using such a sequence guarantees that the approximation error of Equation 2 is not dominated by ∇ ϕ ψM (ϕ). In practice, the target functional ψ(θ) might be quite complex, making gradient estimation challenging. For example, it might involve taking expectations with respect to conditionals of P ϕ (x), or be defined implicitly as a solution to an optimization problem as in [JWZ22a] .One particularly simple and general way to address this challenge is to implement a Monte Carlo estimator of ψ that can be transformed via automatic differentiation into an efficient Monte Carlo estimator for its gradient, a well-understood problem that is beyond the scope of this paper to review. We note that for the very wide class of functionals that can be written as nested expectations, recent work [RCY + 18, SW23, LHSM23] gives formal statements of smoothness assumptions and theoretical results sufficient to obtain the oracle rate ∆ m in terms of numbers of samples, as well as algorithms that can be implemented using automatic differentiation software like PyTorch [PGC + 17]. For example code snippets of functionals, see Appendix B.We conclude by deriving a non-asymptotic error bound for how well Equation 2 approximates Equation 1. For fixed input dimension D and model sizes p, Equation 2 converges to Equation 1 at a O p (1/ √ M ) rate by the Law of Large Numbers. As we are interested in high-dimensional parametric families, we analyze the behavior of our approximation as a function of both input dimension D and model size p. To prove our result, we use standard tools and assumptions from empirical process theory such as the requirement of sub-Gaussian tails [vdV98] . Assumption 3.5. Suppose x ∼ P ϕ (x). There exists a universal constant 0 < C 1 < ∞ such that the normalized score vector x := 1 √ D ∇ ϕ log P ϕ (x) is a sub-Gaussian random vector with parameter C 1 .As E[∥∇ ϕj log P ϕ (x)∥ 2 2 ] = O(D), 1 ≤ j ≤ p, the division by √ D in Assumption 3.5 ensures that the variance of the score does not grow unboundedly as D → ∞. Thus, our assumption that x is sub-Gaussian is mild. Assumption 3.6 below ensures that the functional and score are smooth enough by bounding their gradients. Assumption 3.6. There exist universal constantsC 2 , C 3 < ∞ such that ∥∇ ϕ ψ(ϕ)∥ F < C 2 and ∇ ϕ log P ϕ (x * ) D 2< C 3 for any x * ∈ R D , for any p and D.Unlike our Monte Carlo approximation to the Fisher information matrix, we do not assume a particular type of estimator for ∇ ϕ ψ(ϕ). To prove convergence of MC-EIF to the true EIF, we assume that ∇ ϕ ψM (ϕ) converges to ∇ ϕ ψ(ϕ) at the following rate:Assumption 3.7. Let δ M := ∇ ϕ ψ(ϕ) -∇ ϕ ψM (ϕ) ∈ RL×p denote the approximation error. There exists a universal constant C ψ < ∞ such that for M > C ψ for any ϵ > 0,P ∥δ M ∥ F > p log p + ϵ M < exp(-ϵ), and P ∥∇ ϕ ψM (ϕ)∥ F > C 2 < exp(-ϵ),In Appendix A.1, we prove that Monte Carlo estimators of ∇ ϕ ψ(ϕ) with gradient clipping [ZHSJ20] satisfy Assumption 3.7. Hence, Assumption 3.7 is a mild condition. Under these three assumptions, and the ones in Theorem 3.4, we prove the following result in Appendix A.1, which states that M must scale linearly with p log p to guarantee close pointwise approximation. Theorem 3.8. Suppose Assumption 3.1, Assumption 3.2, Assumption 3.3, Assumption 3.5, Assumption 3.6, and Assumption 3.7 hold. Then, there exists universal constants 0 < C 4 and C 5 < ∞, such that for any ϵ > 0 and M > maxEQUATIONfor x * ∈ R D with probability at least 1 -2 exp(-ϵ), where Σ := cov(x) and λ max (•) denotes the largest eigenvalue of a matrix. 
mc-eif for automated efficient inference: 
 In Theorem 3.8, we proved that MC-EIF is close to the true efficient influence function pointwise. In this Section we; (i) show how MC-EIF can be used to automate the construction of popular efficient estimators, and (ii) prove how many Monte Carlo samples are needed to ensure that key statistical properties hold when MC-EIF is used instead of the true EIF. In doing so, MC-EIF brings conceptual clarity to the practice of constructing efficient estimators, and how these estimators can be implemented using existing differentiable probabilistic programming languages like Pyro [BCJ + 19].All three of the efficient estimator templates we explore in this Section involve some combination of plug-in estimation and EIF-based computations. A key practical benefit of our work is that MC-EIF-based efficient estimators are entirely modular; advances in general-purpose probabilistic inference technology directly translate to advances in efficient estimation under our framework.We start with the simple Von Mises One Step Estimator, which corrects the plug-in estimator in Section 2 by adding the average value of the efficient influence function on a held out dataset. Despite its simplicity, this estimator achieves optimal statistical rates [Ken22] . Our one step estimator using MC-EIF ( φϕ,M (x)) instead of the true efficient influence function (φ ϕ (x)) is provided in Algorithm 1.Input: Target functional ψ, initial estimate of parameters φ, held out datapoints{x n } N n=N/2+1 , Number of Monte Carlo samples M θplug-in ← ψ( φ) {plug-in estimate} C = 2 N N n=N/2+1 φ φ,M (x n ){MC-EIF one step correction} Return: θplug-in + C Theoretical Guarantees. We call the one step estimator that uses the true EIF instead of MC-EIF in Algorithm 1 the analytic one step estimator. Below we prove how many MC samples are needed to ensure our estimator for finite M has the same statistical properties as the analytic one step estimator. Algorithm 2 MC-EIF debiased ML Input: Vector of estimating equations g, initial estimate of parameters φ, held out datapoints{x n } N n=N/2+1 , Number of Monte Carlo samples M f (θ) ← 2 N N n=N/2+1 g(x n , η(p φ), θ) + φ φ,M (x n , θ) {MC-EIF orthogonal moment function} Return: {θ : f (θ) = 0}Construction of Orthogonal Generalized Method of Moment (GMM) Estimators. GMM-based estimators are defined by a nuisance functional η(•) ∈ R J , J ∈ N, and a set of K ∈ N functions {g k (x, η(P ϕ ), θ)} K k=1 , often called estimating equations. These estimating equations are selected so that their roots uniquely identify θ * when the nuisance parameters η(P ϕ ) are estimated correctly:EQUATIONwhereg := (g 1 , • • • , g K ).As an example, g might be the gradient of the log-likelihood function.To make GMM-based estimators less sensitive to incorrect estimation of the nuisance parameters, [CCD + 18, IN22, CNS22] replace g(•) with the orthogonal moment function, constructed using influence functions. In our setting3 , the orthogonal moment function equals the following:EQUATIONwhere φ ϕ (x, θ) is the efficient influence function associated with the functional µ θ (ϕ) = E x∼P ϕ [g(x, η(P ϕ ), θ)] for fixed θ by Equation 2.6 in [IN22] . By Theorem 3.4,EQUATION)Since g is a known by assumption, we can readily use the Monte Carlo methods in [KW14, SHWA15] to automatically approximate ∇ ϕ µ θ (ϕ). We summarize the DML algorithm in Algorithm 2 which replaces Equation 8 with our MC-EIF approximation.Theoretical Guarantees. For general estimating equations, it is difficult to quantity how errors in our MC-EIF approximation to Equation 8 lead to changes in final estimates. When the estimating equations have more structure, however, we obtain a similar result as in Proposition 4.1.Assumption 4.2. g(x n , η(P ϕ ), θ) = m(x n , η(P ϕ ))θ for some vector of known functions m(•). We conclude by writing targeted minimum loss estimation (TMLE) [VDLR06] explicitly in terms of MC-EIF. Unlike the one step estimator or DML, TMLE directly corrects the estimated distribution p φ(x) and then plugs in the corrected distribution into the functional Ψ as the final estimate. To perform this correction it perturbs p φ in the direction of the influence function, searching for the optimal step size by maximizing the perturbed likelihood on the held out dataset. Intuitively, TMLE can be viewed as a form of gradient ascent in function space. We show one step TMLE [VDLR06] in Algorithm 3. The multi-step TMLE version is computed by iterating Algorithm 3 multiple times until ϵ approximately equals 0 [VDLR06] .Algorithm 3 MC-EIF one step TMLEInput: Target functional Ψ, initial estimate of parameters φ, held out datapoints {x n } N n=N/2+1 , Number of Monte Carlo samples M p(ϵ, x) ← (1 + ϵ T φ φ,M (x))p φ(x) {MC-EIF projected ϵ-perturbed density function} ε ← arg max ϵ∈R L :p(ϵ,x)∈M 2 N N n=N/2+1 log p(ϵ, x n ) {Maximum likelihood search over ϵ} Return: Ψ(p(ε, •)) 
experiments: 
 We start by comparing the quality of MC-EIF against other methods for influence function approximation. Then, we show how MC-EIF behaves when; (i) the number of Monte Carlo samples is varied, (ii) the dimensionality of the input is varied, and (iii) the efficient estimator type is varied. Our empirical results ultimately validate our theoretical results in Section 3 and Section 4. Finally, we show how MC-EIF can be used to automate the construction of efficient estimators for new functionals by revisiting a classic problem in optimal portfolio theory. Our MC-EIF implementation is publicly available in the Python package ChiRho. All results shown here are end-to-end reproducible.In [JWZ22a] , the authors target the nonparametric influence function, which is the unique influence function when M = L 2 (P) in Definition 2.1. By contrast, we target the efficient influence function. Thus, for evaluation, we compare how well the empirical Gateaux method from [JWZ22a] approximates the nonparametric influence function and how well our MC-EIF method approximates the efficient influence function on the same data-generating process.To have a ground truth for comparison, we select a simple model and functional where we can analytically compute the nonparametric and efficient influence functions. To this end, we consider the problem of estimating the expected density, Ψ(P) = P(x) 2 dx as in [BR88, CLvdL19] . We further suppose that x ∼ N (µ, σ). We consider two parametric model families: one where µ is unknown but σ = 1, and one where both µ and σ are unknown, which we call M 1 and M 2 respectively. As the nonparametric influence function makes no assumptions on the underlying model family, it remains fixed across M 1 and M 2 and always equals 2(P(X) -Ψ(P)) [CLvdL19] . However, the EIF equals zero for M 1 , as the expected density does not depend on µ. Hence, any plug-in estimate for models in M 1 will result in a correct value of the expected density, and thus no distributional perturbations produce any change. In M 2 , the efficient influence function for the expected density depends on the unknown σ. See Figure 7 in the Appendix for further intuition around the expected density influence functions in parametric (in unknown σ) and non-parametric settings.Figure 1 summarizes how well the empirical Gateaux derivative method approximates the nonparametric influence function and how well our MC-EIF method approximates the EIF at the point Figure 2 : Empirical evidence for convergence theory. Increasing p for the ATE experiments produces MC-EIF approximation errors that closely match and sit below the worst-case error rates given by Theorem 3.8.P ϕ = N (0, 1). We see that MC-EIF is able to approximate the efficient influence function very well (M = 104 samples). By contrast, the empirical Gateaux derivative is highly sensitive to the choice of two kernel smoothing hyperparameters, ϵ and λ. As the true influence function is not known, it is not always clear how to select ϵ and λ. 4 Such numerical instability was already discussed in [CLvdL19] ,where the precision necessary must get exponentially smaller with input dimension, making it infeasible when D ≈ 10.5 MC-EIF, however, has only a single tunable parameter (M ), where larger M unambiguously provides a better approximation. In Theorem 3.8, we provided conditions for this improvement, and Figure 9 of the Appendix corroborates the unsurprising improvement empirically. We further discuss challenges in automating the empirical Gateaux method in Appendix C. We attempted to use the empirical Gateaux derivative as a baseline for other experiments, but were unable to achieve numerically stable solutions for any p > 2 without prohibitively long run-times.Next, we focus on a classic model consisting of a binary treatment, high-dimensional continuous confounders, and Gaussian distributed response; see Appendix E for the precise model formula. We assume that the analyst is interested in estimating the average treatment effect (ATE), where the true ATE is zero but unknown. All influence function computations are relative to an initial point estimate φ, found through maximum a posteriori estimation using 500 training datapoints. Due to the exponential runtime in dimension for the methods in [JWZ22a, CLvdL19] , we focus on comparing MC-EIF with the analytic influence function for ATE below.Sensitivity to Dimensionality. Theorem 3.8 implies that for a fixed number of Monte Carlo samples M , the quality of the approximation degrades with the square root of model dimension p. In Figure 2 , we empirically show how approximation quality degrades as p increases for M = 10 4 fixed. Based on Figure 2 , the empirical results closely match the theoretical behavior predicted by Theorem 3.8. 6 We also show how the computational complexity of MC-EIF scales as p increases in Figure 2 .Sensitivity to Estimator Type. Here, we consider a high-dimensional setup where there are 200 confounders but only 500 training datapoints. We simulate 100 different datasets with this configuration to approximate the sampling distribution of different efficient estimators. In Figure 3 , we see that across estimators, using MC-EIF instead of the true EIF results in minimal downstream error. This is consistent with our theoretical results in Section 4. While MC-EIF is agnostic to the choice of to the oracle estimator. efficient estimator, one may prefer some over others depending on the problem. See Figures 5 and 6 of the Appendix for an example performance comparison between efficient estimators of the ATE.Ability to Handle New Functionals. To illustrate MC-EIF's flexibility, we revisit a classic problem in optimal portfolio theory. Suppose that x ∈ R D is a vector of asset returns. We are interested in estimating the optimal portfolio weights θ * ∈ R D that maximize the expected return while minimizing the variance of the portfolio. Then, the Markowitz optimal portfolio [Mar52] is given by:EQUATION)where λ is the tradeoff between expected return and variance (measure of risk), and Cov(x; P ϕ ) denotes the covariance matrix with respect to P ϕ . Hence, the optimal weights functional Ψ λ (P ϕ ) depends on a high-dimensional nuisance, namely the D × D covariance matrix of returns. The target θ * ϕ,λ = Ψ λ (P ϕ ) is a much lower D-dimensional target parameter. Setting λ = ∞ corresponds to the global minimum variance portfolio [HB91, JM03, ARU20], for which there is (to our knowledge) no efficient estimator in the literature. We show results in Table 1 indicating substantial improvement in a synthetic data experiment; a detailed description of this experiment may be found in Appendix E.2. 
limitations: 
 As discussed in Assumptions 3.1 and 3.2, both the likelihood and the target functional must be differentiable with respect to ϕ. In practice, especially if the model involves latent discrete random variables, some degree of relaxation, marginalization, or reparameterization may be required to ensure differentiability [JGP17] . Recall also that while MC-EIF operates on models with finite parametrizations (Section 3), its capacity to handle high-dimensional nuisance parameters means it can likely apply to, for example, function approximators that recover some of the value proposition offered by non-parametric model components [HSW89] . Additionally, as discussed in Appendix D, infinitedimensional models (like the Gaussian process) can often be reduced to finite ones where MC-EIF can be applied. That said, future work is needed to fully explore the practical and empirical capabilities of MC-EIF in these settings, including how the polynomial complexity of Fisher information matrix inversions plays out in practice. 
conclusion: 
 We have shown both theoretically and empirically that MC-EIF can reliably be used to automate efficient estimation. Our key contributions include MC-EIF's consistency and capability to achieve optimal convergence rates. Empirical evidence shows that MC-EIF performs comparably to traditional estimators using analytic EIFs. Additionally, we illustrate the practical application of MC-EIF in scenarios where the analytic EIF is not known. Given these contributions, there are many exciting areas of future work. For example, one may with to construct more powerful provably efficient estimators on top of MC-EIF (see Appendix D) and explore the growing connection between semiparametric theory and heuristic methods in deep learning [VACB22, BNL + 22, DKSM21, ZDJ + 23]. Additionally, there are many methods that could be used to accelerate the calculation of the Fisher information matrix, which is a computational bottleneck in MC-EIF. Given its foundational role in statistics, various techniques-such as using Kronecker-factored approximations [GM16]-could improve efficiency without sacrificing performance.  
 
 Appendix: EQUATIONBy Assumption 3.7, ∥δ M ∥ F < p+ϵ M with probability greater than 1exp(-ϵ) when M > C ψ . If we can prove that there exists a constant C 4EQUATIONwith probability greater than 1 -2 exp(-ϵ), the claim follows by an application of the union bound. By Theorem 10 in [KK21] , the claim follows if we can prove that there exists a universal constant C 4 such thatEQUATIONwith probability greater than 1exp(-ϵ), where ∥ • ∥ ψ2 denotes the Orlicz sub-Gaussian norm; see Equation 9 in [KK21] for a precise definition of the Orlicz norm of a random vector. With probability greater than 1exp(-ϵ), ∥∇ ϕ ψM (ϕ)∥ F < C 2 by Assumption 3.7. Hence, with probability greater than 1exp(-ϵ),EQUATIONwhere the first inequality follows from Assumption 3.6, the third by [DJS08] and Assumption 3.5, and last by the definition of the spectral norm of a matrix. The result now follows by settingC 4 = C 1 C 2 C 3 .Assumption 3.7 Holds for Monte Carlo Estimators. Here we show if ψM is also approximated with M Monte Carlo samples, then Assumption 3.7 holds. To this end, supposeEQUATIONfor some distribution q(w) and function g ϕ . Such a decomposition exists, for example, when the functional is expressible as a stochastic computation graph [SHWA15] or for reparameterizable densities [KW14] . Suppose further that there exists a universal constant such that ∇ ϕj g ϕ (w m ) ∈ R L is a sub-Gaussian random vector with parameter σ ψ for 1 ≤ j ≤ p. Then,EQUATIONBy Exercise 2.12 in [Wai19] ,max 1≤l≤L,1≤j≤p |[(∇ ϕj g ϕ (w m )] l -[(∇ ϕj ψ(w m )] l | = O p σ ψ log(pL) M (16) Hence, since L is a constant, ∥δ M ∥ F = O p p log(p) M. Thus, the first equation in Assumption 3.7 holds. Under Assumption 3.6, the second equation in Assumption 3.7 trivially holds using gradient clipping with C 2 .A.2 Proof of Theorem 4.1Proof. We want to prove that difference between the analytic one step estimator and Algorithm 1 with finite M decays at a o p 1 √ N rate. Their difference equalsEQUATION)Let ϵ > 0 and supposeM ′ = C 2 4 max(C 5 , 1)(p log(p) + ϵ) max(C 2 1 , 1)N λ 2 max (Σ -1 ) = O(p log pN ),where C 1 is defined in Assumption 3.5, and C 5 and Σ -1 are defined in Theorem 3.8. IfEQUATION)holds, then the proof is complete since M grows faster than M ′ . Let ϵ N = ϵ + log(N/2). By Theorem 3.8 and the union bound,EQUATIONNow,EQUATION)The proof now follows from Equation 19 and Equation 20. Proof. We claim φ ϕ (x, θ) = φ ϕ (x) for all θ. To prove this claim, notice that ∇ ϕ µ θ (ϕ) = ∇ ϕ µ θ ′ (ϕ) for arbitrary θ and θ ′ since µ θ (ϕ) = E x∼P ϕ [m(x n , η(P ϕ ))]θ. Hence, the claim follows from Equation 8. Now,2 N N n=N/2+1 g(x n , η(p φ), θ) + φ φ(x n , θ) = 1 N N n=1 m(x n , η(p φ)) + φ φ(x n ) -θ. (21) Hence, θDML = 2 N N n=N/2+1 m(x n , η(p φ)) + N n=1 φ φ(x n ).Proof of Proposition 4.3. By Lemma A.1, Algorithm 2 uses the same correction term C in Algorithm 1. Hence, the proof of Proposition 4.3 now follows from Proposition 4.1.Remark A.2. By Lemma A.1, the only difference between Algorithm 2 and Algorithm 1 is a different value for the initial estimate of θ * . Specifically, in DML, the initial estimate of θ * is2 N N n=N/2+1 m(x n , η(p φ)), which averages over datapoints drawn from the true distribution. By contrast, Algorithm 1 uses θplug-in , which averages over datapoints simulated from p φ.The implementation of differentiable functional approximations is fairly straightforward when using modern autodifferentation tools. For example, the squared density functional for a mean-zero, univariate normal can be approximated using Monte Carlo as follows:1 N N n=1 N x n , σ 2 ; x n ∼ N 0, σ 2 (22)This can be implemented in pytorch [PGC + 17], and thus automatically differentiated with respect to σ (called "scale" in the code block below) using, for example, torch.autograd.grad. C So, What's Automatic?Figure 4 : We taxonomize the workflow of robust estimation into three stages: the derivation of an (approximate and/or efficient) influence function, the numerical derivation and analysis required for its computation, and the code required to compute it. For the analytic workflow, the derivation of the IF results in Equation 24. This largely involves terms already required by the original plug-in (Equation 23), but still must be implemented on a case-by-case basis in code. For the "Empirical Gateaux" workflow, the first stage requires only the general purpose Equation 25, but demands case-specific numerical considerations and derivations like the one shown in Equation 26. In stark contrast, given a differentiable approximation to the functional of interest, MC-EIF "automates" each stage through use of an end-to-end, general purpose solution.The work required to perform robust estimation can be subdivided into a few key steps. The process begins with a functional of interest, Ψ. With this functional in hand, an analyst must first derive the influence function (or an approximation thereof), and consider any nuances in numerically approximating that quantity. Finally, an engineer must implement that approximation as executable code. Different approaches boast varying levels of "automation" for each step. We claim that in problems where our conditions hold (as outlined in Section 3.2), MC-EIF provides end-to-end automation via a general-purpose solution at each stage. Here, we contrast our approach with both the analytic (see e.g. [Ken16] ) and "Empirical Gateaux" workflows [JWZ22b, CLvdL19] . We will track a workflow's "products" at each of the three stages: first, the derivation of the (efficient and/or approximate) influence function; second, a tractable version of the influence function that properly considers its numerical nuances; third, executable code that computes the influence function. Because our approach uses a general purpose formulation for each of these three stages, we call our approach "automated."We assume the workflow starts having identified a functional of interest and having implemented, in code, a plug-in estimator for it. Throughout, we will use the mean-potential outcome (MPE) functional as our working example7 :EQUATIONThe analytic workflow begins by deriving a closed form influence function -a challenging task even for seasoned experts. This first stage culminates in the following analytic influence function for the MPE [Ken16] .EQUATIONFor general functionals, the derivation resulting in Equation 24 is challenging, even for experts -but given such a derivation, it is often the case that the computation of the quantities composing it can share code and numerical considerations developed to estimate the original, plug-in functional (e.g. Equation 23). Indeed, the influence function of the MPE (Equation 24) involves only terms that an analyst has already considered and implemented for the plug-in (Equation 23). For this reason, we say that in the "analytic" workflow, most of the labor must be allocated to deriving the influence function -tractable, well behaved computation of that influence function tends to involve straightforward extensions of tooling and analysis that already exists for the plug-in.The last stage is the implementation of that tooling in computer code, which will always require some work on a case-by-case basis.The workflow presented by [JWZ22b] and [CLvdL19] significantly reduces the resources required in the first stage -the derivation of the influence function -by providing a general purpose, finite-difference approximation to the influence function (Equation 25). Pϵ,λ , here, represents a perturbation of the estimated distribution P of size ϵ in the direction of a λ-smooth kernel centered at observation O.EQUATIONAt first glance, it seems that computing this term would follow easily given a general purpose framework for the perturbation of P , and then applying the plug-in functional to that perturbed density. Unfortunately, computing Ψ P O ϵ,λ presents a number of numerical challenges in practice. As exhibited in Figure 1 (which echoes figure 1 in the work by [CLvdL19] ), selecting appropriate perturbation parameters ϵ and λ a priori is challenging, and a battle-tested framework for doing so has not yet been developed. Further, in high dimensions (where MC-EIF excels), the required ϵ can be so small as to quickly overrun floating point accuracy on modern computers when even D ≈ 10 [CLvdL19]. Indeed, [JWZ22b] have explicitly left thorough numerical analysis of this approach to further work. In footnote 7, they anecdotally report that quadrature methods were overly sensitive in evaluating perturbed densities in the MPE functional, and instead present a Monte Carlo approach tailored to the task. Unfortunately, neither the numeric considerations or code-implementations of the plug-in estimator easily translate when computing the plug-in with respect to the perturbed data distribution. Below, we show their numeric approximation8 of Equation 25 for the MPE, where observation o comprises (x, a, y), the perturbation kernel K has bandwidth λ, and they use N Monte Carlo samples from a uniform kernel over confounder x.EQUATIONIndeed, just like finite differencing can simplify multivariate calculus, but introduce numeric challenges, the empirical-gateaux approach makes variational calculus easier, but introduces numeric challenges. In sum, we consider the second, "numerical," stage of this workflow to be both labor and expertise intensive, even when the curse of dimensionality does not render it moot.Like in the analytic workflow, the "coding" stage of course requires case-by-case implementations. Moreover, added numerical challenges here introduce significant nuance in implementation that isn't present in the analytic case.In stark contrast, our workflow exploits general solutions in all three stages for the "price" of differentiability of a parametric plug-in estimator. When our general conditions are met (as outlined in Section 3), Equation 1 provides the general purpose solution to the first stage of deriving an (approximate, efficient) influence function, and the second stage is achieved with the Monte Carlo approximation in Equation 2.The third stage is met in software implementing this general purpose solution that operates on functional implementations using one of many auto-differentiation tools now ubiquitous in machine learning (see Appendix B.1 for a simple example). As discussed at the end of Section 3.1, this sometimes requires the ability to exploit methods like the reparameterization trick for the functional of interest. In many cases, modern automatic differentiation software makes this trivial (as shown in Appendix B.1). For some functionals, however, like those involving inner optimizations, this may be more challenging.This general purpose approximation underpins the end-to-end automation of MC-EIF, and is to the best of our knowledge the only such general purpose approximation for efficient influence functions.As a generalization of the gradient operator on ordinary functions, the EIF viewed as an operator on functionals can be shown to have a number of convenient algebraic properties [Ken16, Ken22] , many of which are inherited directly by the MC-EIF estimator. In this section, we speculate on several ways in which these properties could be used to extend the basic MC-EIF framework (and the MC-EIFbased robust estimators in Section 4) to new classes of models and functionals, significantly increasing the range of practical use cases addressable by an implementation of MC-EIF in a differentiable probabilistic programming language like Pyro [BCJ + 19].Multi-argument functionals Many important quantities in statistics and machine learning, like the mutual information I[X; Y ] or KL-divergence KL[P ; Q] are functionals of more than one probability distribution. As shown in [KKP + 15], we can define partial EIFs analogous to partial derivatives for these quantities (which can then be plugged into the efficient estimators of Section 4) by treating all but one argument as part of the functional and computing the ordinary EIF with respect to that argument.Higher-order EIFs Although all of the efficient estimators of Section 4 are derived from the first-order EIF, there are some circumstances where incorporating higher-order EIFs can be shown to be theoretically necessary for achieving certain statistical properties [RLT + 08, BKW23]. Just as ordinary higher-order derivatives are computed by recursively applying a first-order derivative operator to its output, higher-order EIFs can be computed by recursively applying a first-order EIF operator to its own output [vdV14] , a property straightforwardly inherited by MC-EIF.Models with latent variables Thus far, we have assumed that we can exactly simulate from model predictive distributions x ∼ p ϕ and compute log-densities log p ϕ (x), score functions ∇ ϕ log p ϕ (x) and Hessian-vector products. However, our MC-EIF estimator can be extended straightforwardly to models with latent variables and intractable densities and score functions by using a nested Monte Carlo procedure [RCY + 18, SW23] to approximate the prior predictive or posterior predictive distributions and plugging the resulting stochastic estimates into the vanilla MC-EIF framework. We expect our theoretical results to extend to this case provided the approximation error can be made small relative to the Monte Carlo error in estimating the Fisher matrix from a finite set of model Monte Carlo samples.Infinite-dimensional models and targets Semiparametric statistics is by definition fundamentally concerned with models that contain infinite-dimensional (i.e. function-valued) components. There is also intense interest in deriving efficient, doubly robust estimators for infinite-dimensional target functionals like the conditional average treatment effect (CATE) in causal inference. Fortunately, in many of these settings the infinite-dimensional quantities can be reduced to finite ones (and ultimately must be to be representable on a digital computer) to which MC-EIF may be applied in a straightforward way. For example, a Gaussian process is fully characterized by the latent function's values on a finite set of test points; computing the EIF for a functional of the GP reduces to computing the EIF of the finite-dimensional joint distribution on function values at the test points, which is straightforward to estimate with MC-EIF. Similarly, in the case of the CATE, we are ultimately interested in the values of the CATE function on a finite set of test inputs, reducing the problem to ordinary MC-EIF for a finite-dimensional target functional.In Section 5, we consider the following model with confounders c, treatment t, and response y:EQUATIONτ ∼ N (0, 1), (treatment weight) c n ∼ N (0, I D ), (confounders)EQUATIONwhere F ∈ N denotes the number of confounders. In this example, x = (c, t, y) ∈ R D , where D = F + 2 and ϕ = (µ 0 , ξ, π, τ ) ∈ R 2F +2 . To obtain a point estimate in Section 5, we take the maximum a posteriori estimate. In Section 5, we vary the model dimension p by varying F since p = 2F + 2.We assume x is drawn from a multivariate Gaussian distribution with unknown covariance matrix for D = 25 and N = 1000 datapoints. We randomly sample the true covariance matrix using a Lewandowski-Kurowicka-Joe distribution on positive definite matrices. We evaluate MC-EIF and the one step estimator using the relative expected volatility (REV) and the root mean-squared-error (RMSE) between the estimated and the true optimal portfolio weights. Here, the expected volatility is calculated by applying the estimated weights with the actual covariance to the objective in Equation 9. Repeating our experiment using 50 randomly generated datasets, we find that MC-EIF enables substantially improved estimates, as shown in Table 1 .Here, we provide experimental results that provide interesting insight, but do not directly support the key claims of our paper.Figure 5 : Comparison of plug-in estimator and efficient estimators using MC-EIF and analytic EIF for estimating ATE with synthetic data. The true ATE is 0. Closer to zero the better. The distribution is over 100 simulated datasets. Dashed lines represent the estimates using the analytic EIF, and the solid lines represent using MC-EIF (when applicable). Given the high-dimensionality of the problem, the estimation leads to non-zero centering (i.e., some bias remains even after influence function based corrections). Importantly, this is a property of the influence corrected estimators, and is not an artifact introduced by MC-EIF. Instead, we chose our empirical study to demonstrate that MC-EIF produces near-identical results for a diversity of statistical tasks, with a diversity of statistical implications. (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility.In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. Justification: The research presented in this submission does not represent societal impacts except for those shared by all fundamental research. Our work represents general methodological progress for estimating statistical quantities using modern automatic differentiation and probabilistic programming systems. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 11. Safeguards Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? Answer: [NA] Justification: Not applicable Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. 12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected? Answer: [NA] Justification: Not applicable. Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.There are some small variations to this three-step recipe such as using cross-fitting [CCD + 18] instead of a simple equal split of the data; see also Section 4.In Section3.2, we require that ∆m = o m -1 p log p .DML can handle functionals which are not pathwise differentiable. As we only consider parametric model families, however, we can assume without loss of generality that the functional is pathwise differentiable[FS23] and that the range space of the nuisance functional is finite-dimensional. Consequently, we can express DML in terms of efficient influence functions.One suggestion is to visually inspect an epsilon-lambda plot (like the one in Figure1) for a "...possibly curvilinear triangular region nested in the upper left portion of the[plot]. In this region, the finite-difference approximation of the EIF value should be essentially constant ([CLvdL19], §4.3)."Specifically, where the non-parametric dimensionality of the data unit is d = d1 + d2, current theory requires ϵ ≪ λ d 1 , or even ϵ ≪ λ d in some largely avoidable cases ([CLvdL19], §4.1). If functional evaluation is approximate, smaller ϵ requires significantly more compute for the finite difference to dominate Monte Carlo error in estimating the difference between the functional evaluated on the plugin distribution, and the functional evaluated on the ϵ-perturbed distribution. ϵ might even be so small as to overrun floating point accuracy on modern machines -even .1 16 exceeds standard precision recommendations for64-bit floating point values.6 As discussed in Section 3.2, to make the error not grow with p, we would need M ≍ p log p.For simplicity in exposition, but without loss of generality in regards to this description of what is and isn't automated, we follow[JWZ22b] in assuming that outcome and confounders are continuous real numbers, while the treatment is binary.This is their equation 73 in appendix E.3. 


Decision: This paper develops a general Monte-Carlo approach for approximating efficient influence functions. It aims to provide an automated procedure within an existing probabilistic programming workflow. It also provides non-asymptotic error bounds on the quality of the Monte-Carlo approximation and examines three specific estimators: von Mises one step estimator, debiased/double ML, targeted minimum loss estimation. The paper is interesting and fits well with the current developments in the literature. It would be useful to revise the paper for the following aspects: (i) real-data experiments, which the authors indicate that they already carried out; (ii) an emphasis on the advantage in terms of computational speed; (iii) a more serious discussion of limitations (e.g., computation of the Fisher information matrix) as well as the underlying assumptions and requirements.

Reviews: review: paper_summary: paper_summary: The paper proposes Monte Carlo Efficient Influence Functions (MC-EIF), an automated technique for numerically computing efficient influence functions using existing differentiable probabilistic programming systems. MC-EIF simplifies efficient statistical estimation for high-dimensional models, achieving optimal convergence rates and consistency without the need for complex manual derivations. The approach is validated both theoretically and empirically.
strength_weakness: strength_weakness: I am not fully familiar with this line of work, so I am unable to identify a major weakness. However, I have some questions regarding the assumptions  and the effectiveness of the proposed approach on real datasets, which I have added to the Question sections.
questions: questions: Q1. In Assumption 3.5, authors assume that the normalized score vector is sub-Gaussian with a parameter $C_1$. I want to know if this constant scales with respect to the dimension $D$ and model size $p$. If yes, what is the scaling? If no, can you clarify why it does not scale?

Q2. The authors assume that the map $\phi$ to $P_{\phi}$ is continuous. My question is, in order to approximate the Fisher information $\hat{I_{m}}$, do we need to know what $P_{\phi}$ actually is?


Q3. The authors back up their theoretical results with synthetic data experiments. Given the assumptions, it is not clear how applicable the proposed method is to real datasets, and how should the results be interpreted if these assumptions don't hold?"
limitations: limitations: The authors clearly mention the assumptions before stating their theoretical guarantees. Also, the assumptions are explained clearly. I don't see any potential negative societal impact of their work.
paper_summary: paper_summary: The paper establishes a novel method for estimating the efficient influence functions under mild assumptions, called Monte Carlo Efficient Influence Functions(MC-EIF). The method is easy to apply and flexible in many cases, where it can seamlessly equip an existing/popular efficient estimator.
strength_weakness: strength_weakness: When applying new methods in practice, especially in high-dimensional cases, the time cost is also a necessary aspect that needs to be considered, so it would be good to show the time cost for your methods and how it compares to other existing methods.

Some settings in the experiments are not listed clearly, like what is the value of $D$ and $p$ for experiments in Figure 1,2,3 or provided the value of $F$ instead?
questions: questions: It draws my attention that, from Assumption 3.7 and Theorem 3.8, as the model size $p$ becomes larger, the constant of the bound of error will increase as $\sqrt{p \log p}$. So, when dealing with the high-dimensional problem, it may require people to use a really large number of samples ($M$) to get a desirable accuracy and the cost will be tremendous, which could be a potential drawback of this method, what do you think?
limitations: limitations: I am not sure what is the $p$ tested here, but from Figure 2, it seems like $p_{\max}$ is only $1000$, which is not a very large model. People should test on much higher dimensional problems to defend their methods. Moreover, it would be worth looking at more challenging problems, other than well-defined Gaussian problems.
paper_summary: paper_summary: Efficient influence functions (EIFs) for nonparametric estimands are used to construct debiased estimators. Existing methods are primarily estimand-specific and require intricate analytic derivations; existing automated methods don't scale well. This paper proposes general Monte-Carlo estimators of the efficient influence function (MC-EIF) that to be used within a probabilistic programming workflow, shows their convergence to the true efficient influence function, as well as convergence of some estimates based on MC-EIF to the true value.

Disclaimer: I had little familiarity with the use of influence functions for constructing estimators before reading this paper.
strength_weakness: strength_weakness: * I found the paper was harder to understand than was necessary due to unclear notation. In the Problem Statement, P and Q are used interchangably to mean the pdf or the probability measure (e.g. in Definition 2, it's a measure with respect to which the L2 space is defined, and a function in said L2 space).
* The Assumptions (esp 3.1-3.3) are not easily interpretable and are not discussed; as the paper proposes a practical method, it would be made stronger by a discussion of the assumptions.
questions: questions: * The maximum eigenvalue term in Theorem 3.8 may grow in N. Is there an assumptions making sure it does not?
limitations: limitations: * Authors do not discuss if there is value for lifting the assumptions of the method.

Comments: title: Response to rebuttal
Re.
comment: Dear Reviewer hhJt:

Can you please respond to the rebuttal as soon as possible? Your comments will be greatly appreciated. Many thanks,

AC
Thank you very much for your responses! My main concern revolves around the time scale of the methods, which I believe is a key strength. However, this aspect hasn't been sufficiently highlighted in the paper. I hope the author can give this more attention in the revision and also follow through on addressing the larger problem as promised. Overall, I will maintain my score.
Thanks for the response to my questions. A clarifying note on the Continuity assumption, as the authors mentioned, would be very useful. I still believe that the paper would benefit from some real data experiments or results. Overall, I will keep my score for the paper.