Nonparametric Instrumental Variable Regression through Stochastic Approximate Gradients

We thank the reviewer for the insightful observations. Below, we separately address each weakness:

1. The referee is correct in pointing out that estimating conditional densities in high dimensional settings poses nontrivial difficulties, being a limitation which applies in our case, as well as in other approaches to the NPIV problem. In Section 1.1, we start by noting that classical NPIV methods, like those based on sieves or Nadaraya-Watson kernel estimators, exhibit difficulties when dealing with large, high dimensional datasets, while more modern approaches leverage deep learning techniques in an attempt to overcome these issues, but end up susceptible to higher prediction variance under the opposing scenario of few data points. Our goal with this discussion is to point out that, since our formulation is agnostic with respect to how the density ratio is estimated, it could be tailored to specific scenarios and leverage recent advances from other areas. For example, one could use neural networks in large data/dimension scenarios, and kernel methods for situations where there is not a lot of data/the data is low dimensional.

2. We would argue that estimating $\Phi$ is not significantly harder than estimating $\mathcal{P}$, as the latter is an operator acting on an infinite dimensional space. This is evidenced by the very few results discussing convergence rates for estimators of this conditional expectation operator. Furthermore, according to Proposition 3.3., in our loss-agnostic framework, the gradient at $h$ is given by the application of the adjoint operator $\mathcal{P}^*$ at the composition of $\partial_2 \ell$ with $\mathcal{P}[h]$. So, any gradient-based algorithm needs to somehow tackle the estimation of both $\mathcal{P}$ and $\mathcal{P}^*$ as operators. For instance, in [1] they consider sieve estimators using separate basis expansions for $\mathcal{P}$ and $\mathcal{P}^*$. Our stochastic gradient formulation is able to overcome this problem by substituting the outer application of $\mathcal{P}^*$ with a simple multiplication by $\Phi$.

3. To our understanding, in Figure 1, our MSE is lower and statistically better than the standard two-stage method, except for the case in which TSLS is correctly specified. Note, however, that this is expected against any method, as it falls within the realm of parametric estimation. The log-MSE shows significant improvement over the other methods. In particular, the Deep learning based methods (DeepGMM, DeepIV) show larger, more variable log-MSE than the two variants of SAGD-IV. The same is true for Dual IV. The strongest competitor is the KIV method, which demonstrates similar log-MSE for the step and linear function scenarios, lower log-MSE for sine (but with larger variance) and statistically higher log-MSE in the Abs scenario. Therefore, we argue that the numerical experiment indicates that our method is at least comparable with the state-of-the-art methodologies.

   As indicated in our answer to the reviewers first point, we are not claiming that our method is the first to address the high-dimensional setting, we have just pointed out the flexibility of our framework to possibly deal with such scenarios. Given the short reviewing period and the computational resources available to us, unfortunately, we could not run the experiments again in high-dimensional settings.

References

[1] Darolles, S., Fan, Y., Florens, J. P., & Renault, E. (2011). NONPARAMETRIC INSTRUMENTAL REGRESSION. Econometrica, 79(5), 1541–1565. http://www.jstor.org/stable/41237784

We thank the reviewer for providing a careful assessment of our work. In what follows, we address the weaknesses and questions in a point by point fashion.

Weaknesses

1. We are assuming that "consistency" here means convergence of $||\hat{h} - h^*||$ to $0$ as the number of data points grows to infinity. (In [1] the authors use this word to mean convergence of the projected risk to zero, the same type of convergence we provide in Theorem 4.3.).

   Consistency for NPIV regression is, as we mention in Remark 4.5, a very challenging problem. The cited methods [2] and [3], as well as others which prove this kind of result (e.g. [5]), are based on Ridge regression and, hence, can leverage the source condition, ubiquitous in the inverse problems literature, to obtain consistency. In the context of SGD, convergence to the global minimum could be achieved when the function being minimized is strongly convex, which, in our context, translates to restrictions on the operator $\mathcal{P}$, as we point out in Remark 4.5. A similar assumption was also made in [4, Theorem 7, Appendix E.4].

   With respect to sample-complexity, note that $|\mathcal{D}|$ (cf. Assumption 4.2) can be seen as the first stage sample size, while the number $M$ in Algorithm 1 is the second stage sample size. The first two terms on the RHS of Equation (10) quantify the rate of convergence with respect to $M$, while the dependence on $|\mathcal{D}|$ is in the $\zeta$ term. Note that this term is comprised of three estimation errors, whose exact dependence on $|\mathcal{D}|$ can only be known after deciding on specific estimation methods. We did not expand this term further since we chose to be agnostic with respect to the first stage estimators. Nonetheless, if one wants to obtain explicit bounds in $|\mathcal{D}|$, rates for $||\hat{\Phi} - \Phi||$ can be found in [6, Theorem 14.16], while rates for $||\hat{\mathcal{P}} - \mathcal{P}||$ are present in [1, Theorem 2] and [7, Theorem 5]. The other term, $||\hat{r} - r||$, 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$. The book [8] contains several results of this type.

   We agree that this is a discussion which should be in the main text, and we will update Section 4.2 to include it.

2. The novelty of our work does not come necessarily from the proof technique of Theorem 4.3, but rather from the insight that, by tackling NPIV with functional stochastic gradient descent and by carefully exploiting its structure, we are able to obtain an algorithm with improved flexibility that is also capable of addressing other loss functions. This reframing is what allows not only the development of an intuitive and competitive algorithm for the challenging NPIV problem, but also the use of standard tools from convex analysis to theoretically analyze it, which we believe to be an important contribution to the literature and to the NeurIPS community.

   Moreover, while we agree with the referee that the proof technique is standard, we point out that Theorem 4.3 differs from straightforward convex analysis in the sense that the operator associated with the NPIV inverse problem is unknown. This is addressed in Lemma A.3, where we must carefully analyze a term which usually vanishes in simpler instances of statistical inverse problems.

Questions

1. See weakness 2.

2. Thank you for pointing out that more details are needed to parse the algorithm implemention. We will enrich the discussion in Appendix C and the code will be made available. SAGD-IV needs one dataset $\mathcal{D}$ of samples from the triplet $(X, Z, Y)$ to compute $\hat{\Phi}, \hat{\mathcal{P}}$ and $\hat{r}$ (cf. Assumption 4.2 and Remark C.1). Then, it needs another dataset, say, $\mathcal{D}\_Z$, of samples from only $Z$ to conduct the functional gradient descent in Algorithm 1. After one has the estimators $\hat{\Phi}, \hat{\mathcal{P}}$ and $\hat{r}$ in hands, to evaluate $\hat{h}(x)$ one must loop through the samples in $\mathcal{D}\_Z$, computing the stochastic gradient evaluated at $x$ and using it to perform the projected gradient descent step. There are optimizations which can speed up this process, for instance, the term $\partial_{ 2 } \ell \left( \hat{ r } ( z_{ m } ), \hat{ \mathcal{P} } [ \hat{ h }\_{ m - 1 } ] ( z\_{ m } ) \right)$ does not depend on $x$ and, hence, may be left pre-computed. Then, the only effort in computing $\hat{h}(x)$ is in computing $\hat{\Phi}(x, z_m)$ for $1 \leq m \leq M$.

3. We believe that our experiments already fall within a relatively small data regime. Nonetheless, we agree that it is interesting to explore this direction and we ran the experiments with half the sample size. The results can be found in the attached PDF and will be included in the Appendix. We can see that the relative performance of each method largely stays the same. However, as is common with fewer data points, the log-MSE distribution for each method generally exhibited higher variance.

References

[4] Muandet, K., Mehrjou, A., Lee, S. K., & Raj, A. (2020). Dual instrumental variable regression. Advances in Neural Information Processing Systems, 33, 2710-2721.

[5] Darolles, S., Fan, Y., Florens, J. P., & Renault, E. (2011). NONPARAMETRIC INSTRUMENTAL REGRESSION. Econometrica, 79(5), 1541–1565. http://www.jstor.org/stable/41237784

[6] Sugiyama M, Suzuki T, Kanamori T. Density Ratio Estimation in Machine Learning. Cambridge University Press; 2012.

[7] Talwai, P., Shameli, A., & Simchi-Levi, D. (2022, May). Sobolev norm learning rates for conditional mean embeddings. In International conference on artificial intelligence and statistics(pp. 10422-10447). PMLR.

We thank the reviewer for the attentive analysis and relevant suggestions. Below, we separately address the weaknesses and questions.

Weaknesses

Thank you very much for pointing out such an interesting paper. The cited work explicitly uses NN classes as approximations to $L^2$ spaces and conducts SGD on the weights of these neural nets. This is different from our approach, which formulates the SGD updates within the function space $L^2(X)$ and does not assume a parametric form for the stochastic gradients. Another difference is that our work directly addresses the primal problem and Liao et al. (2020) converts it to a saddle-point problem in order to apply a primal-dual algorithm, as you mentioned. Additionally, we do not have any regularization and we do not assume $r$ is known (page 5, paragraph 2, $b$ in their notation). We will include this discussion in the literature review.

Regarding the matter of "risk" versus "projected risk", we wholly agree and will include an appropriate remark in Section 3.

Questions

In [1, Theorem 14.16] the authors provide a convergence rate for our proposed norm of $((\log n) / n)^{1/(2 + \gamma)}$, where $\gamma$ controls the decay of the spectrum of the kernel used to estimate the density ratio.

Concerning the references for the operator regression convergence rate, we thank the reviewer for bringing the second and third to our attention, which will be included in the final version of the paper. We could not find Proposition S5 in the third paper but if the referee can clarify this misunderstanding, we would appreciate it. We will add a paragraph in Subsection 4.2. pointing out these known rates.

References

[1] Sugiyama M, Suzuki T, Kanamori T. Density Ratio Estimation in Machine Learning. Cambridge University Press; 2012.

We thank the reviewer for the constructive comments!

Concerning the weakness, we appreciate the reviewer's input on our exposition and we will improve the discussion on IV and the comparison with current approaches.

Regarding your question, it is possible to use SAGD-IV for other models. However, it depends on finding an appropriate risk functional $\mathcal{R}$, which itself depends on finding an appropriate pointwise loss $\ell$. With the pointwise loss in hands, our method is immediately applicable after computing the derivative $\partial_2 \ell$.

For a general approach, if we assume that $Y = \Lambda(h^*(X), \varepsilon)$ for some nonlinear function $\Lambda$, we have to find a function $F$ such that $\mathbb{E}[ Y \mid Z ] = \mathbb{E}[\Lambda(h^*(X), \varepsilon) \mid Z ] = F( \mathbb{E}[h^*(X) \mid Z] )$. Then, if $Y$ is real-valued, a possible choice for $\ell$ is the $L^2$ loss $\ell(y,y') = \frac{1}{2} (y - F(y'))^2$, and, if $Y$ is discrete, one may use the binary cross entropy loss $\ell(y, y') = BCE(y, F(y'))$, as was done in Section 6 for the binary outcome setting, where $\Lambda(h(X), \varepsilon) = 1 \\{ h(X) + \varepsilon > 0 \\}$.