Nonparametric Instrumental Variable Regression through Stochastic Approximate Gradients

Authors: yuri r fonseca
caio f l peixoto
yuri f saporito
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 17:02:56

Abstract: Instrumental variables (IVs) provide a powerful strategy for identifying causal effects in the presence of unobservable confounders. Within the nonparametric setting (NPIV), recent methods have been based on nonlinear generalizations of Two-Stage Least Squares and on minimax formulations derived from moment conditions or duality. In a novel direction, we show how to formulate a functional stochastic gradient descent algorithm to tackle NPIV regression by directly minimizing the populational risk. We provide theoretical support in the form of bounds on the excess risk, and conduct numerical experiments showcasing our method's superior stability and competitive performance relative to current state-of-the-art alternatives. This algorithm enables flexible estimator choices, such as neural networks or kernel based methods, as well as non-quadratic loss functions, which may be suitable for structural equations beyond the setting of continuous outcomes and additive noise. Finally, we demonstrate this flexibility of our framework by presenting how it naturally addresses the important case of binary outcomes, which has received far less attention by recent developments in the NPIV literature.

Main paper: 
 introduction: 
 Causal inference from observational data presents unique challenges, primarily due to the potential for confounding variables that can affect both outcomes and variables of interest. The unconfoundedness assumption, crucial in this context, posits that all confounding variables are observed and properly accounted for, allowing for an unbiased estimation of causal effects. However, in many real-world scenarios, this assumption is difficult to satisfy. When this is the case, approaches that rely on instrumental variables (IVs) -quantities that are correlated with the variable of interest (relevance condition), do not affect the outcome in any other way (exclusion condition) and are independent of the unobservable confounders -offer a way to identify causal effects even when unobserved confounders exist. As a concrete example, suppose we want to estimate the impact of years of education on earnings. Most likely, there are unobservable factors, such as omitted ability, affecting both the decision to study and income. In this case, changes in compulsory schooling laws could be used as an instrument [2] . performance in high dimensions [35, 25] , Deep Learning based algorithms can suffer from high variance and instability if the amount of available data is small, as was seen in [25, 4] . Furthermore, in many applications of interest, ranging from consumer behavior [30] to epidemiology [23] , the outcome is binary; for a review, see [12] . In these scenarios, a straightforward application of the quadratic loss function would result in a misspecification of the problem, potentially leading to erroneous estimates. The above NPIV methods strongly leverage the additive structure of the problem formulation, and proper extensions would require significant effort, a detailed discussion is presented in Appendix B. Traditionally, a common approach for binary outcomes is using a semiparametric specification based on control function [20] . This allows for nonparametric estimation of the distribution of the unobservables together with identification of the causal function's parameters [1, 32, 14, 11] . To the best of our knowledge, the only work that addresses fully nonparametric IV for binary outcomes is [10] , which proposes an estimation strategy based on Tikhonov regularization.We propose a novel algorithm for NPIV estimation, SAGD-IV, which works by minimizing the populational risk through stochastic gradient descent in a function space. Under mild assumptions, we provide finite sample bounds on the excess risk of our estimator. Empirically, we demonstrate through numerical experiments that SAGD-IV achieves state-of-the-art performance and better stability.We have the freedom to employ a variety of supervised learning algorithms to form the estimator of the stochastic gradient, most notably kernel methods and neural networks. This means our estimator could be tailored to specific scenarios where a particular method is more likely to perform well. Furthermore, our algorithm is naturally able to handle non-quadratic loss functions, which allows us to extend both our algorithm and theoretical guarantees to the binary outcomes case.Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?Answer: [Yes] Justification: We do not overstate the importance of our work and results. The content of the abstract and introduction properly reflect the impacts of our theoretical and experimental results.Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 
problem setup and notation: 
 Following most of the recent literature, we start by presenting the problem setup under the additive noise assumption to demonstrate our methodological contribution. Then, in Section 6, we analyze the important case of binary outcomes, a prototypical example of when this assumption does not hold.Let (Ω, A, P) be the underlying probability space, and let X be a random vector of covariates taking values in X ⊆ R d X . We assume that the response variable Y is generated according toEQUATIONwhere ε ∈ L 2 (Ω, A, P) and satisfies E[ε] = 0. We denote by P X the distribution of the r.v. X and assume that the structural function h ⋆ belongs to L 2 (X , B(X ), P X ) 3 , which we simply denote by L 2 (X). This is a Hilbert space with norm and inner product given by ∥h∥ 2 L 2 (X) = E[h(X) 2 ] and ⟨h, g⟩ L 2 (X) = E[h(X)g(X)]. We assume that E[ε | X] ̸ = 0, that is, some covariates are endogenous. Finally, we assume the existence of a random vector Z, taking values in Z ⊆ R d Z and satisfying 1. E[ε | Z] = 0, the exclusion restriction; 2. X ̸⊥ ⊥ Z, i.e., Z is relevant. This makes Z a valid instrumental variable. We define P Z and L 2 (Z) analogously to P X and L 2 (X). We further consider the mild assumption that X and Z have a joint density denoted by p X,Z . Our goal is to estimate h ⋆ based on i.i.d. samples from the joint distribution of X, Z and Y .As we have listed in the introduction, there are a few different approaches to estimate h ⋆ . We will follow here the original one of [28] , in which we take the expected value of Equation (1) EQUATIONThis motivates us to define the conditional expectation operator P : L 2 (X) → L 2 (Z) given byP[h](z) = E[h(X) | Z = z].This is a bounded linear operator which satisfies ∥P∥ op ≤ 1 (the operator norm) and whose adjoint P * : L 2 (Z) → L 2 (X), that is also a bounded linear operator, is given by P* [g](x) = E[g(Z) | X = x]. Defining r(Z) = E[Y | Z],we can then rewrite Equation (2) as r = P[h ⋆ ].(3) This is a Fredholm integral equation of the first kind [22] and, as such, poses an ill-posed linear inverse problem.In this context, a common assumption [13, 7, 34] made about P is compactness. We also need it here, noting that compact operators with infinite dimensional range provide prototypical examples of ill-posed inverse problems [8] . However, we wish to phrase this assumption in a different, albeit equivalent [8] , form:EQUATIONdenote the ratio of joint over product of marginals densities for the X and Z variables, with the convention that 0/0 = 0. We assume that this kernel has finite L 2 (P X ⊗ P Z ) norm, that is,∥Φ∥ 2 L 2 (P X ⊗P Z ) = X ×Z Φ(x, z) 2 p X (x)p Z (z) dxdz < ∞.This is equivalent to assuming that the conditional expectation operator P : L 2 (X) → L 2 (Z) is Hilbert-Schmidt and, hence, compact.Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes] Justification: Limitations of the algorithm are discussed in Section 4.1.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
the risk and its gradient: 
 Motivated by Equation ( 3), we introduce a pointwise loss function ℓ : R × R → R and define the associated populational risk measure4 R :EQUATIONThe example the reader should keep in mind is the squared loss function ℓ(y, y ′ ) = 1 2 (yy ′ ) 2 , although, as we will see, other examples could be used depending on the particular regression setting. Our goal is to solve the NPIV regression problem by solvinginf h∈H R(h),where H is a closed, convex, bounded subset of L 2 (X) such that h ⋆ ∈ H. We also require 0 ∈ H.For future reference, we state these conditions: Assumption 3.1 (Regularity of H). The set H is a closed, convex, bounded subset of L 2 (X), which contains the origin and satisfies h ⋆ ∈ H.The only part of this assumption which concerns the data generating process is h ⋆ ∈ H, which essentially means that the set H is large enough. 5For H satisfying Assumption 3.1, we let D ≜ diam H < ∞, so that ∥h∥ < D for every h ∈ H. One possible choice for the setEQUATIONwhere A > 0 is a constant. This set is obviously convex and bounded in the L 2 (X) norm. It can be shown that it is also closed, but not necessarily compact, a restriction on the search set imposed in [28] . This can be seen by taking a ∥•∥ ∞ -bounded orthonormal basis for L 2 (X), if one exists. We denote by proj H the orthogonal projection onto H. In case H is given by Equation ( 6), we have the explicit formula6 proj H[h] = (h + ∧ A) -(h -∧ A).We now state all the assumptions needed on the pointwise loss ℓ. We denote by ∂ 2 a partial derivative with respect to the second argument.Assumption 3.2 (Regularity of ℓ).1. The function ℓ : R × R → R is convex and C 2 with respect to its second argument;2. The function ℓ has Lipschitz first derivative with respect to the second argument, i.e., there exists L ≥ 0 such that, for all y, y ′ , u, u ′ ∈ R we have|∂ 2 ℓ(y, y ′ ) -∂ 2 ℓ(u, u ′ )| ≤ L(|y -u| + |y ′ -u ′ |).Some useful facts which follow immediately from these assumptions are presented in Appendix A.Here we divert from the previous literature and proceed to compute a functional stochastic gradient for the risk R. A similar idea has been considered by [15] in the context of statistical inverse problems, however, their setup assumes that the operator posing the inverse problem is known, which greatly simplifies the analysis and cannot be directly applied to the NPIV problem. We start by providing an analytical formula for ∇R(h):Proposition 3.3. The risk R is Fréchet differentiable and it's gradient satisfiesEQUATIONwhereP * : L 2 (Z) → L 2 (X)is the adjoint of the operator P.Both P and P * are unknown operators which commonly appear in NPIV estimation [8, 13] , with Nadaraya-Watson kernels [26, 39] being a classical option for estimating them. However, the fact that they are nested in Equation ( 7) is undesired, and may impose extra computational costs, especially considering that ∂ 2 ℓ is nonlinear in the second argument for non-quadratic losses. We overcome this difficulty by leveraging the following characterization of the gradient:Corollary 3.4. The gradient of the populational risk satisfiesEQUATIONwhere Φ is defined as in Equation (4).The benefit of our approach is that Equation (8) enables a more computationally efficient way to estimate the gradient of the risk functional, as we explain next. From Equation ( 8), for a given x ∈ X , the random variable Φ(x, Z)∂ 2 ℓ(r(Z), P[h](Z)) is an unbiased stochastic estimator of ∇R(h)(x). Our stochastic approximate gradient is then built using estimators Φ, r and P of Φ, r and P respectively. With this notation, given a sample Z, we considerEQUATIONWe have then substituted the estimation of the operator P * by the simpler problem of estimating the ratio of densities Φ and computing its product with ∂ 2 ℓ( r(Z), P[h](Z)). Moreover, density ratio estimation is an area of active research within the machine learning community, which makes developments in this direction immediately translatable into benefits to our method.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
algorithm: implementation and theory: 
 In Algorithm 1 we present Stochastic Approximate Gradient Descent IV (SAGD-IV) 7 , a method for estimating h ⋆ using the approximation given by Equation ( 9).Input: Samples (z m ) M m=1 . Estimators Φ, r and P. Sequence of learning rates (α m ) M m=1 . Initial guess h 0 ∈ H.Output: h for 1 ≤ m ≤ M do Set u m = Φ(•, z m )∂ 2 ℓ r(z m ), P[ h m-1 ](z m ) Set h m = proj H h m-1 -α m u m end for Set h = 1 M M m=1 h m Remark 4.1.We note the fact that the internal loop of Algorithm 1 only needs samples from the instrumental variable Z to unfold.For our theoretical analysis, all we require of the estimators of Φ, r and P is the following: Assumption 4.2 (Properties of Φ, r and P).1. Φ, r and P are computed using a dataset D, composed of N samples of the triplet (X, Z, Y ), which are independent from the Z samples used in Algorithm 1.2. r ∈ L 2 (Z) a.s.;3. P : L 2 (X) → L 2 (Z) is a bounded linear operator a.s. 4. ∥ Φ∥ ∞ ≜ ess sup (x,z)∈X ×Z | Φ(x, z)| < ∞. This implies, in particular, ∥ Φ∥ L 2 (P X ⊗P Z ) < ∞.Algorithm 1 is purposefully agnostic to how the estimators P, Φ and r are obtained, in order to provide the user with sufficient modeling flexibility. In what follows, we will explain the options we considered for accomplishing these tasks in the experiments section. A more detailed description is given in Appendix C.Estimating Φ is a interesting problem in itself as it is a ratio of densities, and it has drawn significant attention from the Machine Learning community. Here, we consider kernel and neural network methods, regarded as the state of the art. On the other hand, not many options are available to compute P, which means estimating all possible regression of functions of X over Z. Here, we chose to use kernel mean embeddings, as was developed [34] . Finally, estimating r is the simplest procedure in our method as it is a regression of Y over Z. Similarly to Φ, we also consider kernel and neural network methods.Based on these options, we decided to specify two variants of SAGD-IV built upon how we compute Φ and r: Kernel SAGD-IV uses RKHS methods for both estimators, while Deep SAGD-IV employs neural network estimators in these two tasks. The method for computing P is kernel mean embeddings in both variants.Under a big data scenario, the dependence of our algorithm on kernel methods to estimate the P operator could be a limitation. However, we are able to avoid this issue for the estimation of Φ and r. Nonetheless, since our algorithm is agnostic to the specific estimator of P, we could take advantage of any recent developments in the literature of operator-estimation problems.Since we are directly optimizing the projected populational risk measure, we are able to provide guarantees for R( h) in mean with respect to the training data z 1:M = {z 1 , . . . , z m }. Our main result is the following: EQUATIONwhereξ = 3 2 ∥ Φ∥ 2 ∞ C 2 0 + L 2 ∥ r∥ 2 L 2 (Z) + L 2 D 2 ∥ P∥ 2 op , ζ = ∥Φ -Φ∥ 2 L 2 (P X ⊗P Z ) + ∥r -r∥ 2 L 2 (Z) + ∥P -P∥ 2 op , τ = 2D max 3(C 2 0 + L 2 E[Y 2 ] + L 2 D 2 ), 2L 2 ∥ Φ∥ 2 ∞ , 2L 2 D 2 ∥ Φ∥ 2 ∞ .It is productive to analyze the RHS of the bound in Equation ( 10) more carefully. If we choose the sequence (α m ) to satisfy the usual assumptionsM α M → ∞ and 1 M M m=1 α m → 0as M → ∞, then the first two terms in the sum vanish as M grows. The last term appears due to the fact that we do not know Φ, r nor P, but it explicitly quantifies how the estimation errors of Φ, r and P come together to determine the quality of the final estimator. Furthermore, it converges to zero if the chosen estimation methods are consistent. 8 Remark 4.4 (Y versus r(Z)). Since r = P[h ⋆ ], one should compute the risk as in Equation ( 5), comparing P[h](Z) with r(Z). However, most works on NPIV [28, 34, 25] compareP[h](Z) directly with Y instead of r(Z).In Appendix E we discuss this difference and its relationship with the estimator r. Remark 4.5 (Consistency for h * ). Since the NPIV regression is an ill-posed inverse problem, consistency of our algorithm is a very challenging and interesting research question, and it would require additional assumptions. For instance, under strong convexity of the populational risk measure, it is well known that consistency is obtained if the excess risk converges to zero. It can be shown that, for the quadratic loss, strong convexity of the populational risk is equivalent to P having a continuous inverse. 
numerical experiments: 
 Here, we compare the performance of the two variants of SAGD-IV explained in Section 4.1, Kernel SAGD-IV and Deep SAGD-IV, with that of five other notable algorithms for estimating h ⋆ in the context of continuous responses: KIV, DualIV, DeepGMM, DeepIV and 2SLS. Hyperparameters and other implementation details are discussed in Appendix C.To study the performance of our estimator in a continuous response setting, we used the data generating process from [6] , which we recall below:EQUATIONZ = (Z 1 , Z 2 ) ∼ Uniform([-3, 3] 2 ),ϵ ∼ N (0, 1), γ, δ ∼ N (0, 0.1).Thus, the confounding variable is ϵ and within this specification, four options for h ⋆ are considered:step : h ⋆ (x) = 1{x > 0}, abs : h ⋆ (x) = |x|, linear : h ⋆ (x) = x, sin : h ⋆ (x) = sin(x).8 If one desires to substitute the ζ term for estimator-specific bounds, rates for ∥Φ -Φ∥ L 2 (P X ⊗P Z ) can be found in [36, Theorem 14.16] , while bounds for ∥P -P∥op are present in [34, Theorem 2] , [38, Theorem 5] and [33, Proposition S5] . The term ∥r -r∥ L 2 (Z) is the risk for a simple real-valued regression problem, and various rates are available depending on the chosen regressor and the degree of smoothness assumed on r [16] .We adjusted each estimator using 3000 random variable samples (see Remark C.1), and tested them on 1000 samples of X. For each method and response function, we evaluated predictions over 20 realizations of the data. Log mean squared error (MSE) box plots and plots of each method's estimator for a randomly chosen realization of the data are displayed in Figure 1 . The first thing we can notice from Figure 1 is that TSLS dominates the linear scenario, as expected, but behaves poorly in all others, showing the importance of nonlinear approaches.On a second note, both variants of SAGD-IV performed competitively in all scenarios, with the KIV method being a strong competitor in step and sin. However, its performance was significantly worse for the abs response function. Additionally, it is worth noticing the competitive MSEs of the deep learning variant of SADG-IV, which performed consistently better than the other two neural-network-based algorithms.Furthermore, an important observation is that, within kernel methods, Kernel SAGD-IV had the best stability/performance tradeoff, having the smallest variance across scenarios. In contrast, DualIV exhibited the highest variance. This can be explained by considering DualIV's hyperparameter selection algorithm [25, Section 4] . In order to perform cross-validation to choose the best regularization parameters, the algorithm requires another (different) regularization parameter as an input, which must be chosen without any guidance. The authors state that it was set to a "small constant". This can explain the poor performance, given the importance of regularization parameters for ridge regression algorithms. We note our method does not suffer from this problem. 
binary response: 
 Here we outline how our methodology can be applied to a scenario where the quadratic loss is not the natural option: that of binary response models. We show that, by simply modifying the loss ℓ in the risk definition and, consequently, in the algorithm, we are able to attack this problem as it is currently formulated in the literature. This is in clear contrast with recent methods for NPIV estimation, which, when applicable, would require significant effort. We discuss the limitations of current methods and possible extensions in Appendix B.In the binary response model we consider, the only change we need to make to the data generating process is the following:Y = 1 {Y ⋆ > 0} ,whereY ⋆ = h ⋆ (X) + ε.That is, we observe a signal which indicates if the response variable is greater than the threshold zero. Although simple, this type of model captures important applications as, for instance, in economics/consumer behavior, where Y * is the utility the consumer gets from buying a product and we only observe Y , the buy/no buy decision. The function h ⋆ captures the impact of the variable X on the consumers' utility. As usual in the literature, h ⋆ is interpretable relative to the scale parameter of ε. Here, X, Z and ε satisfy the same assumptions as before.In order to use one of the available NPIV regression methods in the literature, we must first identify a suitable risk functional, of the form shown in Equation (5) .Define η = h ⋆ (X) -E[h ⋆ (X) | Z] + ε, so that E[η | Z] = 0 and Y = 1 {E[h ⋆ (X) | Z] + η > 0} .Denoting by F η|Z the conditional distribution function of η given Z, we findE[Y | Z] = P [E[h ⋆ (X) | Z] + η > 0 | Z] = 1 -F η|Z (-E[h ⋆ (X) | Z]).Assume that the distribution of η given Z is always symmetric around zero, so that 1 -F η|Z (-y) = F η|Z (y). We then have r(Z) = F η|Z (P[h ⋆ ](Z)). Notice that any notion of risk for a given h ∈ L 2 (X) would have to compare r(Z) and P[h](Z) through F η|Z , which is an unknown function.Hence, without making further assumptions, it is not clear how to setup the risk in a useful way, which hinders current NPIV methods from tackling the problem.As far as we know, the only assumption in the literature which allows one to be fully nonparametric in h ⋆ within the binary outcomes setting was first introduced in [10] . It posits that η is independent of Z, and has known distribution functionF . Since η = Y ⋆ -E[Y ⋆ | Z],it is already uncorrelated with any function of Z. The assumption extends this to independence, and makes available the resulting distribution function.Under this assumption, we conclude that r(Z) = F (P[h ⋆ ](Z)) and, by the negative log-likelihood of the Bernoulli distribution, a natural candidate for ℓ isℓ(y, y ′ ) = BCE(y, F (y ′ )),where BCE is the binary cross entropy function:EQUATIONWhereas it is obvious that the quadratic loss satisfies Assumption 3.2, it is not clear whether the same is true for the loss ℓ(y, y ′ ) = BCE(y, F (y ′ )). Nevertheless, if F (x) = σ(x) ≜ 1/(1 + exp(-x/β)), the c.d.f. of a logistic distribution with scale parameter β, one can verify that Assumption 3.2 holds.We mimicked the continuous response DGP making the necessary modifications:Y = 1{E[h ⋆ (X) | Z] + η > 0}, X = Z 1 + η + γ, Z = (Z 1 , Z 2 ) ∼ Uniform([-3, 3] 2 ), η ∼ Logistic(0, β) γ ∼ N (0, 0.1).Here, β is a scale parameter which was set to √ 0.1. For response functions, we analyzed the sin and linear cases, since we must compute E[h ⋆ (X) | Z] analytically in order to generate the data. This computation is easy for the selected scenarios, and we haveE[X | Z = z] = z 1 , E[sin(X) | Z = z] = βπ • e -0.1 2 sinh(βπ) sin(z 1 )As with the continuous response setting, we present log MSE results and samples from the resulting estimators in Figure 2 .We can immediately see that our kernel based estimator performed remarkably well for both response functions. The binary setup is significantly more challenging, as we are reducing the continuous response Y to a simple binary signal. Kernel SAGD-IV obtained results which are, on average, on par with the continuous response scenario and Deep SAGD-IV had more difficulty uncovering the true structural function. 
conclusion: 
 In this work, we proposed SAGD-IV, a flexible framework for performing nonparametric instrumental variable regression. It is based on a novel way to tackle the NPIV problem, in which we formulate the minimization of the populational risk as functional stochastic gradient descent. We are then able to naturally consider non-quadratic loss function and incorporate various regression methods under one formulation. Under mild assumptions, we provided bounds on the excess risk of our estimator. Furthermore, we empirically demonstrated superior stability and state-of-the-art MSEs for continuous outcomes, as well as a promising performance for the binary regression setup. 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines?Answer: [Yes]Justification: We have conformed with the NeurIPS Code of Ethics.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [NA] Justification: Our work discusses only theoretical and computational aspects of an estimation problem. We do not foresee any societal impact.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). 
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: The paper poses no such risks. 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: [Yes] Justification: All methods are properly cited. 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 authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. 
 
 Appendix: Proposition A.1. Assume that ℓ satisfies Assumption 3.2. Then:1. Setting C 0 = |∂ 2 ℓ(0, 0)|, we have |∂ 2 ℓ(y, y ′ )| ≤ C 0 + L(|y| + |y ′ |), for all y, y ′ ∈ R; 2. The map f → ∂ 2 ℓ(r(•), f (•)) from L 2 (Z) to L 2 (Z)is well-defined and L-Lipschitz.3. The second derivative with respect to the second argument is bounded:∂ 2 2 ℓ(y, y ′ ) ≤ L for all y, y ′ ∈ R. Proof. 1. Write ∂ 2 ℓ(y, y ′ ) = ∂ 2 ℓ(y, y ′ ) -∂ 2 ℓ(0, 0) + ∂ 2 ℓ(0, 0)and apply the triangle inequality as well as Assumption 3.2 Item 2.2. From the previous item we know this map is well-defined. If f and g belong to L 2 (Z), we have∥∂ 2 ℓ(r(•), f (•)) -∂ 2 ℓ(r(•), g(•))∥ 2 L 2 (Z) = E |∂ 2 ℓ(r(Z), f (Z)) -∂ 2 ℓ(r(Z), g(Z))| 2 ≤ L 2 E |f (Z) -g(Z)| 2 = L 2 ∥f -g∥ 2 L 2 (Z) .3. It follows directly from Assumption 3.2 Item 2.Proof of Proposition 3.3. We start by computing the directional derivative of R at h in the direction f , denoted by DR[h](f ):DR[h](f ) = lim δ→0 1 δ [R(h + δf ) -R(f )] = lim δ→0 1 δ E [ℓ(r(Z), P[h + δf ](Z)) -ℓ(r(Z), P[h](Z))] = lim δ→0 1 δ E [ℓ(r(Z), P[h](Z) + δP[f ](Z)) -ℓ(r(Z), P[h](Z))] = lim δ→0 1 δ E δ∂ 2 ℓ(r(Z), P[h](Z)) • P[f ](Z) + δ 2 2 ∂ 2 2 ℓ(r(Z), P[h + θf ](Z)) • P[f ](Z) 2 = E [∂ 2 ℓ(r(Z), P[h](Z)) • P[f ](Z)] + lim δ→0 E δ 2 ∂ 2 2 ℓ(r(Z), P[h + θf ](Z)) • P[f ](Z) 2 = E [∂ 2 ℓ(r(Z), P[h](Z)) • P[f ](Z)] ,where θ ∈ (0, δ) is due to Taylor's formula. The last step is then due to Proposition A.1 Item 3. We can in fact expand the calculation a bit more, as follows:DR[h](f ) = E [∂ 2 ℓ(r(Z), P[h](Z)) • P[f ](Z)] = ⟨∂ 2 ℓ(r(•), P[h](•)), P[f ]⟩ L 2 (Z) = ⟨P * [∂ 2 ℓ(r(Z), P[h](•))], f ⟩ L 2 (X) .This shows that R is Gâteaux-differentiable, with Gâteaux derivative at h given by DR[h] = P * [∂ 2 ℓ(r(•), P[h](•))].By Proposition A.1 Item 2, we have that h → DR[h] is a continuous mapping from L 2 (X) to L 2 (X), which implies that R is also Fréchet-differentiable, and both derivatives coincide [29] . Therefore,EQUATIONProof of Corollary 3.4. From Proposition 3.3 we know that∇R(h)(x) = P * [∂ 2 ℓ(r(•), P[h](•))](x) = E[∂ 2 ℓ(r(Z), P[h](Z)) | X = x] = Z p Z|X (z | x)∂ 2 ℓ(r(z), P[h](z)) dz = Z p X,Z (x, z) p X (x) ∂ 2 ℓ(r(z), P[h](z)) dz = Z p Z (z) • p X,Z (x, z) p X (x)p Z (z) ∂ 2 ℓ(r(z), P[h](z)) dz = E[Φ(x, Z)∂ 2 ℓ(r(Z), P[h](Z))],where we defineΦ(x, z) = p X,Z (x, z)/(p X (x)p Z (z)) if p Z (z) > 0 and Φ(x, z) = 0 if p Z (z) = 0.Notice we may assume that p X (x) > 0, since the set of x ∈ X such that this happens has probability 1, and ∇R(h), being an element of L 2 (X), is only defined almost surely.Before presenting the proof of Theorem 4.3, we need to prove two auxiliary lemmas. To lighten the notation, the symbols ∥•∥ and ⟨•, •⟩, when written without a subscript to specify which space they refer to, will act as the norm and inner product, respectively, of L 2 (X).Lemma A.2. In the procedure of Algorithm 1 we have u m ∈ L 2 (X) for all 1 ≤ m ≤ M and, furthermore,E z 1:M [∥u m ∥ 2 ] ≤ ρ( Φ, r, P), where ρ( Φ, r, P) = 3∥ Φ∥ 2 ∞ C 2 0 + L 2 ∥ r∥ 2 L 2 (Z) + L 2 D 2 ∥ P∥ 2 op .Proof. By Assumption 4.2 we have:∥u m ∥ 2 L 2 (X) = ∥ Φ(•, z m )∂ 2 ℓ r(z m ), P[ h m-1 ](z m ) ∥ 2 L 2 (X) = E X Φ(X, z m )∂ 2 ℓ r(z m ), P[ h m-1 ](z m ) 2 ≤ ∂ 2 ℓ( r(z m ), P[ h m-1 ](z m )) 2 ∥ Φ∥ 2 ∞ < ∞.Hence, u m ∈ L 2 (X) for all m. This computation and Proposition A.1 Item 1 then implyE z 1:M ∥u m ∥ 2 ≤ 3∥ Φ∥ 2 ∞ C 2 0 + L 2 ∥ r∥ 2 L 2 (Z) + ∥ P[ h m-1 ]∥ 2 L 2 (Z) ≤ 3∥ Φ∥ 2 ∞ C 2 0 + L 2 ∥ r∥ 2 L 2 (Z) + ∥ P∥ 2 op ∥ h m-1 ∥ 2 ≤ 3∥ Φ∥ 2 ∞ C 2 0 + L 2 ∥ r∥ 2 L 2 (Z) + D 2 ∥ P∥ 2 op = 3∥ Φ∥ 2 ∞ C 2 0 + L 2 ∥ r∥ 2 L 2 (Z) + L 2 D 2 ∥ P∥ 2 op ≜ ρ( Φ, r, P).Lemma A.3. In the procedure of Algorithm 1 we have∥E zm [∇R( h m-1 ) -u m ]∥ ≤ κ( Φ) ∥Φ -Φ∥ 2 L 2 (P X ⊗P Z ) + ∥r -r∥ 2 L 2 (Z) + ∥P -P∥ 2 op 1 2 , where κ 2 ( Φ) ≜ 2 max 3(C 2 0 + L 2 E[Y 2 ] + L 2 D 2 ), 2L 2 ∥ Φ∥ 2 ∞ , 2L 2 D 2 ∥ Φ∥ 2 ∞ .Proof. To ease the notation, we defineΨ m (Z) ≜ ∂ 2 ℓ(r(Z), P[ h m-1 ](Z)), Ψ m (Z) ≜ ∂ 2 ℓ( r(Z), P[ h m-1 ](Z)).Let's expand the definition of ∥•∥:EQUATION, Now we add and subtract Φ(X, Z)Ψ m (Z), so thatEQUATION,where ∥Φ∥ 2 L 2 (P X ⊗P Z ) = X ×Z Φ(x, z) 2 p(x)p(z) dxdzis the norm with respect to the independent coupling of the distributions of X and Z. By Proposition A.1.1 we have∥Ψ m ∥ 2 L 2 (Z) = E Z ∂ 2 ℓ(r(Z), P[ h m-1 ](Z)) 2 ≤ E Z C 0 + L |r(Z)| + P[ h m-1 ](Z) 2 ≤ 3 C 2 0 + L 2 ∥r∥ 2 L 2 (Z) + L 2 ∥P[ h m-1 ]∥ 2 L 2 (Z) ≤ 3 C 2 0 + L 2 E[Y 2 ] + L 2 D 2 . It is also clear that, by Assumption 4.2, ∥ Φ∥ 2 L 2 (P X ⊗P Z ) ≤ ∥ Φ∥ 2 ∞ .Finally, by Assumption 3.2.2 we also have∥Ψ m -Ψ m ∥ 2 L 2 (Z) = E Z ∂ 2 ℓ(r(Z), P[ h m-1 ](Z)) -∂ 2 ℓ( r(Z), P[ h m-1 ](Z)) 2 ≤ 2L 2 ∥r -r∥ 2 L 2 (Z) + ∥(P -P)[ h m-1 ]∥ 2 L 2 (Z) ≤ 2L 2 ∥r -r∥ 2 L 2 (Z) + D 2 ∥P -P∥ 2 op .To combine all terms, we first defineκ 2 ( Φ) ≜ 2 max 3(C 2 0 + L 2 E[Y 2 ] + L 2 D 2 ), 2L 2 ∥ Φ∥ 2 ∞ , 2L 2 D 2 ∥ Φ∥ 2 ∞ .Then, it's easy to see that∥E zm [∇R( h m-1 ) -u m ]∥ ≤ κ( Φ) ∥Φ -Φ∥ 2 L 2 (P X ⊗P Z ) + ∥r -r∥ 2 L 2 (Z) + ∥P -P∥ 2 op 1 2 ,as we wanted to show.We now have everything needed to show the Proof of Theorem 4.3. We start by checking that R is convex inH: if h, g ∈ H and λ ∈ [0, 1], then R(λh + (1 -λ)g) = E[ℓ(r(Z), P[λh + (1 -λ)g](Z))] = E[ℓ(r(Z), λP[h](Z) + (1 -λ)P[g](Z))] ≤ λE[ℓ(r(Z), P[h](Z))] + (1 -λ)E[ℓ(r(Z), P[g](Z))] = λR(h) + (1 -λ)R(g).By Assumption 3.1, h ⋆ ∈ H9 . By the Algorithm 1 procedure, we have1 2 ∥ h m -h ⋆ ∥ 2 = 1 2 proj H h m-1 -α m u m -h ⋆ 2 ≤ 1 2 ∥ h m-1 -α m u m -h ⋆ ∥ 2 = 1 2 ∥ h m-1 -h ⋆ ∥ 2 -α m ⟨u m , h m-1 -h ⋆ ⟩ + α 2 m 2 ∥u m ∥ 2 .After adding and subtracting α m ⟨∇R( h m-1 ), h m-1h ⋆ ⟩, we are left with1 2 ∥ h m-1 -h ⋆ ∥ 2 -α m ⟨u m -∇R( h m-1 ), h m-1 -h ⋆ ⟩+ α 2 m 2 ∥u m ∥ 2 -α m ⟨∇R( h m-1 ), h m-1 -h ⋆ ⟩.Applying the first order convexity inequality on the last term give us, in total,1 2 ∥ h m -h ⋆ ∥ 2 ≤ 1 2 ∥ h m-1 -h ⋆ ∥ 2 -α m ⟨u m -∇R( h m-1 ), h m-1 -h ⋆ ⟩ + α 2 m 2 ∥u m ∥ 2 -α m (R( h m-1 ) -R(h ⋆ )).Hence, making this substitution and rearranging terms, we getR( h m-1 ) -R(h ⋆ ) ≤ 1 2α m ∥ h m-1 -h ⋆ ∥ 2 -∥ h m -h ⋆ ∥ 2 + α m 2 ∥u m ∥ 2 -⟨u m -∇R( h m-1 ), h m-1 -h ⋆ ⟩.Finally, summing over 1 ≤ m ≤ M leads toEQUATIONThe next step is to take the average of both sides with respect to z 1:M , taking advantage of the independence between z 1:M and D, the data used to compute Φ, r and P. Each summation in the RHS is then bounded separately.The first summation admits a deterministic bound. By assumption, the diameter D of H is finite. HenceEQUATIONThe second summation can be bounded with the aid of Lemma A.2:EQUATIONFinally, the third summation can be bounded using Lemma A. 3 . Let E z-m denote the expectation with respect to z 1 , . . . , z m-1 , z m+1 , . . . , z M and notice thatE z 1:M ⟨∇R( h m-1 ) -u m , h m-1 -h ⋆ ⟩ = E z-m E zm ⟨∇R( h m-1 ) -u m , h m-1 -h ⋆ ⟩ = E z-m ⟨E zm ∇R( h m-1 ) -u m , h m-1 -h ⋆ ⟩ ≤ E z-m ∥E zm ∇R( h m-1 ) -u m ∥∥ h m-1 -h ⋆ ∥ ≤ DE z-m ∥E zm ∇R( h m-1 ) -u m ∥ .Then, applying Lemma A.3 and setting τ ≜ Dκ we getEQUATION)All that is left to do is to apply equations ( 14), ( 15), ( 16) and ( 17) along with the inequality which defines convexity. Let h ≜ 1 M M m=1 h m-1 and ξ ≜ ρ/2. Then:E z 1:M R( h) -R(h ⋆ ) ≤ 1 M M m=1 E z 1:M R( h m ) -R(h ⋆ ) ≤ D 2 2M α M + ξ( Φ, r, P) 1 M M m=1 α m + τ ( Φ) ∥Φ -Φ∥ 2 L 2 (P X ⊗P Z ) + ∥r -r∥ 2 L 2 (Z) + ∥P -P∥ 2 op 1 2 .In this section, we provide some theoretical comparisons with other notable methods for NPIV estimation.In our formulation, we substitute the task for estimating P * by that of estimating the ratio of densities Φ. In this way, we bypass the need to evaluate an operator estimator P * in every iteration, replacing this by evaluating Φ once and multiplying it by the derivative of the loss function (cf. Algorithm 1). Aside from the computational advantage this poses, it is interesting to investigate the differences which arise in theoretical bounds when comparing with a method which ends up using an estimator of P * .In KIV, the authors propose a first stage procedure which estimates E (P, in our notation) and, in the second stage, they use the adjoint of this operator through the kernel mean embedding µ(z) = E * ϕ(z) (ϕ is the feature map, see the notation in [34] ). To estimate E * , they use the adjoint of their estimator E n λ of E, where n is the sample size and λ is a regularization parameter. In our Theorem 4.3, provided the learning rate has been appropriately chosen, the excess risk is upper bounded by ∥Φ -Φ∥ + ∥r -r∥ + ∥P -P∥. In fact, ∥r -r∥ is just the risk for a simple regression problem, and so we will focus our analysis on the other two terms. In KIV, the excess risk is bounded using [37, Theorem 2] , according to which it is smaller than the sum of five terms. We reproduce this bound here for a more detailed analysis: Theorem B.1 (Proposition 32 of [34] ). The excess risk of KIV's stage 2 estimator H m ξ can be bounded by five terms:EQUATIONwhereS -1 = ∥ √ T • ( T + ξ) -1 ( g -g)∥ 2 HΩ , S 0 = ∥ √ T • ( T + ξ) -1 • (T -T)H m ξ ∥ 2 HΩ , S 1 = ∥ √ T • (T + ξ) -1 (g -TH ρ )∥ 2 HΩ , S 2 = ∥ √ T • (T + ξ) -1 • (T -T)(H ξ -H ρ )∥ 2 HΩ , A(ξ) = ∥ √ T (H ξ -H ρ )∥ 2 HΩ .In our notation, E corresponds to R, H m ξ corresponds to h and H ρ corresponds to h ⋆ . The space H Ω is a RKHS which contains the solution H ρ . It is also important to recall the definition of the terms in bold, since they depend on µ n λ (z) = (E n λ ) * ϕ(z):T = 1 m m i=1 T µ(zi) , g = 1 m m i=1 Ω µ(zi) ỹi T = 1 m m i=1 T µ n λ (zi) , g = 1 m m i=1 Ω µ n λ (zi) ỹi .Where zi and ỹi are second-stage samples.We refer the reader to [34] for a detailed description of the other mathematical objects. Here, we will argue thatEQUATIONThe other terms in Equation ( 18) are controlled by optimally choosing how the regularizing parameter ξ goes to zero. This relates to the first two terms of our bound in Equation ( 10) that goes to zero as the SGD iteration increases, as long as we properly choose the learning rate. In [34, Proposition 35] , these two terms are bounded as follows:S -1 ≤ ∥ √ T • ( T + ξ) -1 ∥ 2 L(HΩ) ∥ g -g∥ 2 HΩ , S 0 ≤ ∥ √ T • ( T + ξ) -1 ∥ 2 L(HΩ) ∥ T -T∥ 2 L(HΩ) ∥H m ξ ∥ 2 HΩ .Then, [34, Proposition 37] refers to [37, Supplements 7.1.1 and 7.1.2] to bound ∥ g -g∥ 2 and ∥ T -T∥ 2 . In [37] , one can see that these bounds amount to10 ∥ g -g∥ 2 ≤ C • 1 m m i=1 ∥µ n λ (z i ) -µ(z i )∥ 2 H X |ỹ i | 2 ∥ T -T∥ 2 ≤ D • 1 m m i=1 ∥µ n λ (z i ) -µ(z i )∥ 2 H X ,where C and D are constants. Then, KIV uses the bound on ∥µ n λ (z i )µ(z i )∥ 2 H X provided in their Corollary 1, which is based on the inequality11 ∥µ n λ (z i ) -µ(z i )∥ 2 H X ≤ ∥E n λ -E ρ ∥ 2 HΓ ∥ϕ(z i )∥ 2 H ZHence, with high probability and ignoring constants, we have12 ∥ g -g∥ ≲ ∥E n λ -E ρ ∥ ∥ T -T∥ ≲ ∥E n λ -E ρ ∥.This would be an intermediate result in [34, Proposition 37] .Therefore, the estimation error of E n λ when compared to E impacts the risk of KIV's estimator through two terms. These are analogous to our two terms ∥Φ -Φ∥ and ∥P -P∥, since ∥Φ -Φ∥ L 2 (P X ⊗P Z ) is the Hilbert-Schmidt norm of the difference between P * and the integral operator induced by Φ [8, Theorem 2.42].In Section 6, we show that, under the binary outcomes setting, the correct formulation of the risk isR(h) = E[BCE(r(Z), F (P[h](Z)))].In KIV's notation, this would become13 E(H) = E[BCE(Y, F (Hµ(Z)))]Substituting the populational mean by the empirical one, taking into account that we must estimate µ and applying regularization, we get (cf. [34, Section 4.2] )E m ξ (H) = 1 m m i=1 BCE(ỹ i , F (Hµ n λ (z i ))) + ξ∥H∥ 2 HΩ .Using the fact that H Ω and H X are isometrically isomorphic [34, Section 4.2], we can formulate the risk in an equivalent way, now defined on H X :E m ξ (h) = 1 m m i=1 BCE(ỹ i , F ((E n λ h)(z i )) + ξ∥h∥ 2 H X .Then, the KIV estimator is defined ash m ξ = arg min h∈H X E m ξ (h).For the quadratic loss, KIV finds the solution in closed form, since this is a ridge regression problem. However, for this new loss, no obvious closed form solution exists. Applying the Representer Theorem, we find that the solution h m ξ should satisfyh m ξ = n i=1 α i ψ(x i ),for some α ∈ R n , where ψ(x i ) = K(x i , •) is the reproducing kernel of H X . Then, one is left with the following optimization problem:arg min α∈R n 1 m m i=1 BCE(ỹ i , F (α ′ w(z i )) + ξα ′ K XX α Since (E n λ h m ξ )(z) = α ′ w(z) (see the definition of w in [34, Appendix A.5.1]). An approximate solution to this n-dimensional optimization problem must be found with numerical methods. Keep in mind that this problem must be solved several times in order to find a good value of the regularization parameter ξ via cross-validation, which is crucial for the performance of regularized regression methods. In contrast, our method does not require additional work after defining the modified risk.In DualIV, the risk minimization problem is reformulated into a saddle point optimization problem, which we write in their notation:min f ∈F max u∈U Ψ(f, u),whereΨ(f, u) = E XY Z [f (X)u(Y, Z)] -E Y Z [ℓ ⋆ Y (u(Y, Z))].Here, F is the space which contains the solution to the NPIV estimation problem, U is the space of the "dual function" and ℓ ⋆ y is the Fenchel dual of the function ℓ y = ℓ(y, •). In the binary outcomes setting, the loss is given by ℓ(y, y ′ ) = BCE(y, F (y ′ )). If F is the logistic function, i.e. F (y ′ ) = exp(y ′ )/(1 + exp(y ′ )), one may verify thatℓ ⋆ y (u) = -H(u + y), if 0 ≤ u + y ≤ 1, +∞, otherwise,whereH(p) = -[p log p + (1 -p) log(1 -p)]is the entropy. The empirical minimax objective then becomesΨ(f, u) = 1 n n i=1 f (x i )u(y i , z i ) + 1 n n i=1 H(y i + u(y i , z i )),with the constrain that 0 ≤ y i + u(y i , z i ) ≤ 1 for every i. Again, in contrast with the case of quadratic loss, when one may choose F and U to be RKHS and obtain a closed form ridge-regularized solution, this formulation is not as tractable.In GMM based methods for continuous outcomes, the main idea is that sinceE[ε | Z] = 0 and Y = h ⋆ (X) + ε, we have E[Y -h ⋆ (X) | Z] = 0 and, hence, E[f (Z)(Y -h ⋆ (X))] = 0 for every f ∈ L 2 (Z).If we choose a large number of functions f 1 , . . . , f m , we may search for h ⋆ by finding a function which minimizesEQUATION)This is possible because the error term ε can be written as Yh ⋆ (X). Under the structural equation Y = 1{h ⋆ (X) + ε > 0}, such a characterization of the error term is not available, which already presents a difficulty. Nonetheless, as we show in Section 6, under this new equation we haveE[Y | Z] = F (P[h ⋆ ](Z)). Therefore, moment conditions could be obtained through E[(Y -F (P[h ⋆ ](Z)))f (Z)] = 0 for every f ∈ L 2 (Z).However, this would create the need to introduce an estimator P of P. This defeats one of the main advantages of turning to GMM, which is to avoid the first stage of two-stage methods, since errors in the first stage can have a large impact on the resulting estimator.In this section, we give practical guidelines for implementing SAGD-IV, as well further details regarding the baseline methods below:• KIV: Since the original implementation was only publicly available in Matlab, to facilitate reproducibility of our results we re-implemented this method in Python. This was done following the guidelines presented in the paper. • DualIV: Similarly to KIV, we rewrote the method in Python using the closed form solutions and the hyperparameter selection procedure present in the paper. • DeepGMM: We used the publicly available Python implementation 14 .• DeepIV: We used the implementation available as part of the EconML package [5] .• 2SLS: We include the classic Two Stage Least Squares procedure to show the benefits of following a nonparametric approach.All experiments were conducted on a Apple Sillicon M1 Pro CPU. No GPUs were used. Remark C.1 (Experiments sample size). All four baseline methods we are considering need a single dataset (X, Z, Y ) to fit their estimator (hyperpameter validation included), whereas our method needs one dataset of (X, Z, Y ) triplets to compute Φ, r and r, and one dataset of Z samples to conduct the SAGD loop. To make fair comparisons, we trained each method using the same total amount of random variable samples. For example, DeepGMM might use 800 triplets (X, Z, Y ) for training and 200 triplets for validation, whereas SAGD-IV may use 600 triplets for fitting the preliminary estimators and then 1200 samples of Z for the SAGD loop, both totaling 3000 random variable samples used to adjust the model.C.1.1 On computing P, Φ and r Ratio of Densities. Density ratio estimation is a well studied problem in Machine Learning, with [36] being a good resource for acquainting oneself with relevant algorithms. The approaches we follow here consist on obtaining an estimator Φ that minimizes the (empirical version of the) quantity1 2 ∥ Φ -Φ∥ 2 L 2 (P X ⊗P Z ). We consider two options for Φ: a Reproducing Kernel Hilbert Space (RKHS) estimator obtained through the Unconstrained Least Squares Importance Fitting (uLSIF) framework described in [36] and also a neural network trained to minimize the same objective over its parameters. By changing the kernel or the activation function used, we can make sure this estimator adheres to Assumption 4.2 Item 4.Conditional expectation operator. This is the most complex object we must estimate, since it encompasses all possible regressions of functions of X over Z. Not many options exist for this task, with a popular one [13, 8] being the use of Nadaraya-Watson kernels [26, 39] . We, however, chose to employ the estimation method used as a first stage in KIV [34] , which transforms the problem into vector-valued RKHS ridge regression, due to its straightforward implementation.Conditional mean of Y given Z. This, in contrast, is the simplest unknown in our formulation, being the regression function of Y over Z. We again consider two options for this task: the first is to employ the method previously discussed in order to estimate the conditional expectation operator of Y given Z, and then apply it to the identity function. The second option is to train a general purpose regressor, such as a neural network.Since the final estimator in Algorithm 1 is an average over the whole path of estimates produced by the SAGD loop, we can improve upon its quality by discarding the first K estimates h 1 , . . . , h K , where K must be chosen independently from M . In this way, it is as if we started our algorithm from h K+1 , after a warm up period.Regarding sample splitting, our theoretical results do not specify an optimal ratio between the M SAGD iterations and the N samples of (X, Z, Y ) used to compute Φ, r and P. We empirically determined that M = 2N is a good guideline to follow, and it is the ratio we employed in Section 5. Of course, in real data scenarios this ratio might be constrained due to limited data availability.• The Algorithm 1 learning rate was set to be α m = 1/ √ M for 1 ≤ m ≤ M ;• The warm up time K was set to 100;• The set H was chosen as in Equation ( 6), with A set to 10. During our experiments, we found that this parameter, aside from needing to be large enough for us to have h ⋆ ∈ H, did not have any noticeable influence on the resulting estimator.• For Kernel SAGD-IV:-We used Gaussian kernels, with lenghtscale parameter determined through the median heuristic; -The regularization parameters for the three kernel ridge regressions were chosen through cross validation. Since these parameters are essential for actual learning to take place, we selected the values to be tested through an iterative procedure, described in Algorithm 2.• For Deep SAGD-IV:-For both networks in the roles of Φ and r, we used two hidden dense layers of sizes 64 and 32. In the continuous response setting the activation function used was a ReLU, while for the binary response one we used the sigmoid activation. For Φ, we used dropout with rate 0.01. -Both networks were trained with a batch size of 512, using the Adam [21] optimizer with learning rate of 0.01, accross 1.5 • 10 5 /N epochs. We adopted L 2 regularization in the two estimators, with regularization parameters of 3 • 10 -3 for r and 5 • 10 -3 for Φ. We also implemented early stopping in the two training procedures.• We remind the reader that both SAGD-IV variants use kernel methods to compute P.The modifications we made with respect to the original paper were the following:• Instead of using 50% of the samples for training and 50% for validation, we adopted a train / validation split of 80 / 20.• We changed the batch size from 1024 to 256. This was done because we supplied the algorithm with 1000 training samples instead of the 2000 used in [6] , and hence a batch size of 1024 would not produce stochastic gradients.We followed all the guidelines presented in the original paper for low dimensional settings.We employed a 50/50 split between (X, Z, Y ) and ( X, Z, Ỹ ) observations (see [34] , Appendix A.5.2). Regularization parameter candidates were also chosen using the selection method in Algorithm 2. We followed the hyperparameter selection guidelines given in Section 4 of [25] . This means we used a 50/50 split for train/validation samples and chose the regularization parameters from a grid search, minimizing the squared norm of the dual function u.This section presents a version of the experiments in Section 5.1 with half the sample size, providing a comparison between SAGD-IV and other methods within a smaller data availability regime. Figure 3 contains the corresponding log-MSE results.Contrasting Figures 1 and 3 , we can see that the relative performances for each method mostly stay the same. Some methods, such as Deep-IV in the abs scenario, exhibited an error distribution with larger variance, but this is to be expected from the reduction in sample size, which makes any given dataset realization a worse depiction of the true data distribution. These results support the conclusion that our method's performance is competitive with other state-of-the-art approaches in more data availability conditions. Computing the risk using r(Z) or Y are equivalent under the quadratic loss, however, one might think that considering Y instead of r(Z) in our formulation would lead to an algorithm that does not need to estimate r(Z). Below we show that this is not the case, since the gradients of the risk measure under both formulations are the same. This implies that we would still need to estimate r for the algorithm.If instead of Equation ( 5), we consider the risk defined asR(h) = E[ℓ(Y, P[h](Z))],then proceeding in exactly the same way as in the proof of Proposition 3.3, we would arrive atD R[h](f ) = E[∂ 2 ℓ(Y, P[h](Z)) • P[f ](Z)]As it is, this expression is not an inner product in L 2 (Z) since Y is not measurable with respect to Z. However, if ∂ 2 ℓ is linear with respect to its first argument, which is true for the quadratic and binary cross entropy losses, we haveD R[h](f ) = E[∂ 2 ℓ(Y, P[h](Z)) • P[f ](Z)] = E [E[∂ 2 ℓ(Y, P[h](Z)) • P[f ](Z) | Z]] = E [∂ 2 ℓ(E[Y | Z], P[h](Z)) • P[f ](Z)] = E [∂ 2 ℓ(r(Z), P[h](Z)) • P[f ](Z)] ,which is the same expression obtained before. Hence, the gradient of R still depends on r(Z) and the resulting method still needs to estimate this quantity.Nonetheless, one could still use the sample y i as an estimate of r(z i ). However, this is not the best choice, because the term ∥ r -r∥NeurIPS Paper ChecklistReproducing Kernel Hilbert SpaceWe denote by B(X ) the Borel σ-algebra in X .We note that a more fitting name for this object would be projected risk measure, since we are projecting h onto the instrument space before applying the loss function. Concerning this matter, we chose to follow the terminology used in the current NPIV literature, which refers to this projected risk simply as "risk".This part could be generalized to (h ⋆ +ker P)∩H ̸ = ∅, so that there exists h ∈ H such that R(h) = R(h ⋆ ), without prejudice to any of the theoretical results.Here we use the notation h + = max{h, 0}, h -= (-h) + and a ∧ b = min{a, b}.The nomenclature Stochastic Approximate Gradient Descent (SAGD) first appeared in another paper addressing an unrelated topic,[31], where Langevin dynamics are employed to obtain approximate gradients when samples of the underlying random variable are not available.We want to remind the reader that it is enough to assume that there exists h ∈ (h ⋆ + ker P) ∩ H. By the definition of h, we have R( h) = R(h ⋆ ).In[37] there is a 2h in the exponent instead of 2, which is translated to 2ι in KIV. However, this is only a matter of assuming Holder continuity of Ω with exponent ι. If one assumes Lipschitz continuity, the exponent becomes 2.Eρ = E.We use the symbol ≲ to mean that the LHS is smaller or equal than a constant times the RHS, where the constant does not depend on data.Since BCE is linear in the first argument, it does not matter if we use Y or r(Z) to define the risk.https://github.com/CausalML/DeepGMM.gitL 2 (Z) which appears in the RHS of Equation(10) would becomeE[Var[Y | Z]]. This is constant with respect to the size of the auxiliary dataset D (cf. Assumption 4.2) and could be improved by taking r to be a better regressor of Y over Z. 


Decision: The majority of reviewers tend to agree that the paper provides a solid contribution to the research area of Nonparametric Instrumental Variable (NPIV) Regression. It proposes to solve NPIS using a functional stochastic gradient descent algorithm, which is novel. The theoretical results are also very thorough. I believe that people working in this area will appreciate this work.

Reviews: review: paper_summary: paper_summary: This paper introduces a novel IV method, called the SAGD-IV, which is more efficient and stable in the NPIV regression. Two different variants of the SAGD-IV are given, and a range of comparisons between these variants with the existing methods are given. Moreover, the performance of SAGD is not only shown in a continuous response case but also is extended to a discrete case.
strength_weakness: strength_weakness: There are a lot of mathematical equations and deductions, like theories, assumptions etc, which is great and indicates that you have a pretty solid work, but it is also great to spend a bit more room for storytelling, which will bring people straight to the main point and get your amazing work.
questions: questions: I wonder if the error term is not additive and but is a non-linear function of $X$ and $\epsilon$, would the SAGD-IV methods still work?
limitations: limitations: N/A
paper_summary: paper_summary: This paper considers the standard problem of non-parametric instrumental variable estimation (NPIV) and proposes a new approach of functional stochastic gradient descent to solve it (SAGD–IV), where the gradient estimator can be implemented and adapted using certain machine learning or deep learning techniques. Theoretical guarantee of the finite-time convergence of the algorithm is provided under suitable assumptions. Simulated experiments show the effectiveness of the proposed method to some extent.
strength_weakness: strength_weakness: 1. The consistency or sample-complexity of the proposed method is unknown. This type of finite sample results are presented in recent previous methods on NPIV including but not limited to [1, 2, 3]. 

2. The technique and idea behind the proof of the convergence of SAGD-IV is relatively standard and straightforward given that the risk function $\mathcal{R}$ is convex in the functional space.


**References:**

[1] Singh, R., Sahani, M., & Gretton, A. (2019). Kernel instrumental variable regression. Advances in Neural Information Processing Systems, 32.

[2] Bennett, A., Kallus, N., Mao, X., Newey, W., Syrgkanis, V., & Uehara, M. (2023, July). Minimax Instrumental Variable Regression and $ L_2 $ Convergence Guarantees without Identification or Closedness. In The Thirty Sixth Annual Conference on Learning Theory (pp. 2291-2318). PMLR.

[3] Li, Z., Lan, H., Syrgkanis, V., Wang, M., & Uehara, M. (2024). Regularized DeepIV with Model Selection. arXiv preprint arXiv:2403.04236.
questions: questions: 1. Could you please highlight the technical difficulty and the novelty in proving the main theoretical results of the convergence of SAGD-IV?

2. It is not clearly stated how to exactly perform the functional gradient descent step as proposed in Algorithm 1 in the practical experiments. Could the authors elaborate on this a bit more? 

3. What is the performance of the proposed algorithm in the small data regime when compared with previous methods?
limitations: limitations: Please see the above weakness section and the question section.
paper_summary: paper_summary: The authors propose an estimator for nonparametric instrumental variable regression (NPIV), with an extension to the binary outcome case. They prove that the excess risk of the projected estimator is controlled by rates of the density ratio, a regression, and the conditional expectation operator. Extensive comparative simulations are conducted. Careful analysis in the appendix also compares this approach with others.
strength_weakness: strength_weakness: Algorithm 1 of ''Provably Efficient Neural Estimation of Structural Equation Model: An Adversarial Approach’’ by Liao et al. (2020) is another stochastic approximate gradient approach for NPIV. I believe the difference in the algorithms amounts to whether the primal or dual formulation is optimized with stochastic gradients. A clarification would be welcome.

It would be good to put into words that the excess risk being analyzed is for objects projected onto the instrument space, which is common in this literature but perhaps less so in the NeurIPS community.
questions: questions: Are there known rates for the density ratio given in the first term of zeta? Please provide at least one reference where a rate in this norm, or a sufficient norm, is given.

Similarly, it would be good to point to rates for the third term of zeta, such as Theorem 2 of ''Kernel Instrumental Variable Regression'' by Singh et al. (2019), Theorem 5 of ''Sobolev Norm Learning Rates for Conditional Mean Embeddings'' by Talwai et al. (2022), and Proposition S5 of ''Kernel Methods for Causal Functions: Dose, Heterogeneous and Incremental Response Curves'' by Singh et al. (2024).
limitations: limitations: NA
paper_summary: paper_summary: This paper studies the nonparametric instrumental variable (NPIV) regression problem. The authors present a new algorithm, SAGD-IV, which utilizes stochastic gradient descent in a function space to directly minimize the populational risk. The gradient can be computed via computing the estimating the conditional density $\Phi$ and the conditional expectation operator $\mathcal{P}$. This approach is distinguished by its flexibility, accommodating a variety of supervised learning algorithms, such as neural networks and kernel methods, and supporting both continuous and binary outcome scenarios. Finally, the authors also proved the consistency of the NPIV estimator under regularity assumptions.
strength_weakness: strength_weakness: I have several concerns of the paper: 

1 The author claims that existing methods in estimating NPIV solutions face the problem of cannot apply to large, high dimensional datasets. Typically in the NPIV context, such a problem arises due to the estimation of the condtional density function $\Phi$. However, in the paper, computing the gradient also requires us to compute the conditonal density fucntion, so I do not see how such limitation on high dimensional dataset can be avoided.

2. The novelty of the algorithm relies on performing the sgd in the second stage after we estimate $\Phi$ and $\mathcal{P}$. I believe in standard two stage estimation methods, once we obtain $\mathcal{P}$ (no estimationf of $\Phi$ required), the conditional expectation operator, we could also apply standard gradient descent to obtain the final estimator. I do not see how the proposed algorithm is significantly different from the standard ones. However, the proposed algorithm have the additional requirement to compute the conditional density function which is significantly harder than computing $\mathcal{P}$.

3. The conducted experiments is not convicing as even in the low dimensional regime, the obtained mse is not significantly better than the standard two stage estimation method. And the experiments do not contain the high dimensional setting. So it is not clear how the method works in such settings.  A nolvety the authors seem to claim.
questions: questions: See weakness.
limitations: limitations: NA

Comments: title: End of rebuttal period
Response to rebuttal
comment: As the rebuttal period comes to an end and all reviewers have replied, we would like to thank everyone for taking the time to read and respond to our rebuttals. We are specially grateful for the comments which regard our work as new, interesting, and relevant to the NeurIPS community.

The feedback given by the reviewers has motivated a series of improvements to our paper, including:

* A more detailed discussion on the decay rate for the $\zeta$ term in Theorem 4.3, with the addition of new references for bounding the estimation errors of $\Phi$, $\mathcal{P}$ and $r$;

* More details on the algorithm's implementation and another experiment in a smaller data size setting;

* An improvement on the discussion of IV regression and related works in the introduction.
Thank you very much for your detailed response to all my questions and concerns! Also I appreciate the additional experiments conducted in the more small data regimes. Given the theoretical contributions and the potentially insights this work could bring to the community, we would update my score to 5.
Thanks for those replies! I will maintain my score.
Thanks for the clarification, we will add the reference in the final version of the paper.
Thank you for these replies.

Proposition S5 in the third paper can be found in its Online Supplement on the Biometrika website.
Thanks for the reply. I will maintain my score.