Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
59,567 characters · 12 sections · 58 citation commands
Fast and Adaptive Rates for Regularized DeepIV
\twocolumn[ \icmltitle{Fast and Adaptive Rates for Regularized DeepIV}
\icmlsetsymbol{equal}{*}
\icmlaffiliation{yyy}{Princeton University, NJ, USA} \icmlaffiliation{comp}{Genentech, CA, USA} \icmlaffiliation{sch}{Stanford University, CA, USA}
\icmlcorrespondingauthor{Zihao Li}{[email removed]}
\icmlkeywords{Machine Learning, ICML}
\vskip 0.3in ]
\printAffiliationsAndNotice{\icmlEqualContribution} \newcount\Comments \Comments=1 \definecolor{darkgreen}{rgb}{0,0.5,0} \definecolor{darkred}{rgb}{0.7,0,0} \definecolor{teal}{rgb}{0.3,0.8,0.8} \definecolor{orange}{rgb}{1.0,0.5,0.0} \definecolor{purple}{rgb}{0.8,0.0,0.8} \usetikzlibrary{positioning}
\definecolor{Gray}{gray}{0.9}
Instrumental variable (IV) estimation is an important problem in various fields, such as causal inference angrist1995identification, newey2003instrumental,deaner2018proxy,cui2020semiparametric,kallus2021causal,kallus2022causal, missing data problems miao2018confounding,wang2014instrumental, dynamic discrete choice models kalouptsidi2021linear and reinforcement learning liao2021instrumental,uehara2022future,uehara2022provably,shi2022minimax,wang2021provably,yu2022strategic.
In this paper, we focus on nonparametric IV (NPIV) regression newey2003instrumental. NPIV concerns three random variables $X \in \mathbb{R}^d$ (covariate), $Y \in \mathbb{R}$ (outcome variable), and $Z \in \mathbb{R}^d$ (instrumental variables). We are interested in finding a solution $h_0$ of the following conditional moment equation dikkala2020minimax, chernozhukov2019inference: $$ \mathbb{E}[Y - h(X)|Z] = 0. $$ This is equivalently written as ${\mathcal{T}} f=r_0$ where ${\mathcal{T}}: L_2(X)\ni f(X) \mapsto \mathbb{E}[f(X)|Z]\in L_2(Z)$ and $r_0(Z) = \mathbb{E}[Y|Z]$ by denoting $L_2(X), L_2(Z)$ to be the $L_2$ space defined on $X$ and $Z$ with respect to the underlying distribution. Both the operator ${\mathcal{T}}$ and $\mathbb{E}[Y|Z]$ remain unknown. Hence, we aim to solve ${\mathcal{T}} f=r_0$ by harnessing an identically independent distributed (i.i.d.) dataset $\{X_i, Y_i, Z_i\}_{i\in[n]}$.
There has been a surge in interest in NPIV regressions that try to integrate general function approximation such as deep neural networks beyond classical nonparametric models hartford2017deep,singh2019kernel,xu2021deep,zhang2023instrumental,dikkala2020minimax,bennett2020variational,bennett2023minimax,bennett2023source,kallus2022causal,singh2020kernel. Despite these extensive efforts, existing approaches encounter several challenges. The first challenge is the ill-posedness of the inverse problem. Many existing works liao2020provably,newey2003instrumental,florens2011identification,kato2021learning require that the NPIV solution $h_0$ is unique, and further impose quantitative bounds on measures of ill-posedness. However, it is known that the uniqueness assumption is easily violated in practical scenarios, such as weak IV andrews2005inference, andrews2019weak or proximal causal inference kallus2021causal. The second challenge involves the reliance on minimax optimization oracles in many methods bennett2023minimax,dikkala2020minimax,liao2020provably,bennett2023source,zhang2023instrumental, which results in minimax non-convex non-concave optimization when invoking deep neural networks. However, currently, such an optimization can be notoriously unstable and may fail to converge lin2020near,jin2020local,lin2020gradient,diakonikolas2021efficient,razaviyayn2020nonconvex. Instead, our approach seeks to address this challenge by proposing a computationally efficient estimator that relies on standard supervised learning oracles rather than minimax oracles. The third challenge is the absence of clear procedures for model selection in existing works xu2021deep,zhang2023instrumental,cui2020semiparametric,hartford2017deep. This issue is problematic, because model selection, including techniques like cross-validation, has played a pivotal role in the practical success of machine learning algorithms bartlett2002model,gold2003model,guyon2010model,cawley2010over,raschka2018model,emmert2019evaluation,mcallester2003pac. Model selection becomes essential particularly in scenarios where the true NPIV solution $h_0$ lies outside the chosen function classes optimized by the algorithm, which has been seldom explored in prior works.
To address aforementioned challenges, we propose a two-stage method, which we refer to as the Regularized DeepIV (RDIV). This approach consists of two steps. First, we learn the operator ${\mathcal{T}}$ by maximum likelihood estimation (MLE). Secondly, we obtain an estimator for $h_0$ by solving a loss incorporating the learned ${\mathcal{T}}$ and Tikhonov regularization ito2014inverse to handle scenarios where solutions of the conditional moment constraint are nonunique. While our method can be viewed as a regularized variant of the DeepIV method of hartford2017deep with a non-parametric MLE first-stage, no prior theoretical convergence guarantees exist for the DeepIV method. We show that our estimators can converge to the least norm IV solution (even if solutions are nonunique) and derive its $L_2$ error rate guarantee based on critical radius. Subsequently, we introduce model selection procedures for our estimators. Particularly, we provide theoretical guarantees for model selection via out-of-sample validation approaches, and show an oracle result in our context. Finally, we further illustrate that our method can be easily generalized to an iterative estimator that more effectively leverages the well-posedness of $h_0$.
Our contribution is to propose the first estimator for NPIV that (a) operates in the absence of the uniqueness assumption, (b) does not rely on the minimax computational oracle, and (c) allows for model selection. Subsequently, we demonstrate that our estimator can be extended to an iterative estimator, which achieves a state-of-the-art convergence rate in terms of $L_2$ error analogous to bennett2023source, while bennett2023source requires a minimax computational oracle and does not permit us to perform model selection. Therefore, our estimator can be seen as an estimator with a strong theoretical guarantee due to the property (a) while it is practical due to properties (b) and (c). Notably, none of the existing works can enjoy such a guarantee, as shown in Table (ref).
For a function $f: \mathcal{X} \times \mathcal{Y} \times \mathcal{Z} \rightarrow \mathbb{R}$, we denote its population expectation by $\mathbb{E}[f(X,Y,Z)]$. We denote the empirical mean of $f$ by $\mathbb{E}_n[f(X,Y,Z)]:= \frac{1}{n}\sum_{i=1}^n f(X_i, Y_i, Z_i)$. We denote the set of all probability distributions defined on set $\Omega$ by $\Delta(\Omega)$. We denote the $L_p$ norm of $f$ by $\|f\|_p := \mathbb{E}[|f|^p]^{1/p}$. Throughout the paper, whenever we use a generic norm of a function $\|f\|$, we will be referring to the $L_2$-norm. For two density function $p(x)$ and $q(x)$, we denote their Hellinger distance by $H(p(\cdot)\mid q(\cdot)) = \int_\mathcal{X} (\sqrt{p(x)} - \sqrt{q(x)})^2 d\mu(x)$. For a functional operator ${\mathcal{T}}: L_2(X) \rightarrow L_2(Z)$, we denote the range space of ${\mathcal{T}}$ by $\mathcal{R}({\mathcal{T}} )$, i.e., $\mathcal{R}({\mathcal{T}}) = \{{\mathcal{T}} h: h \in L_2(X)\}.$ Moreover, we use ${\mathcal{T}}^*: L_2(Z) \rightarrow L_2(X)$ to denote the adjoint operator of ${\mathcal{T}}$, i.e.,$\langle g, {\mathcal{T}} h\rangle_{L_2(Z)} = \langle {\mathcal{T}}^* g, h\rangle_{L_2(X)}$ for any $h \in L_2(X), g \in L_2(Z)$, where $\langle\cdot,\cdot\rangle_{L_2(X)}$ and $\langle\cdot,\cdot\rangle_{L_2(Z)}$ are inner products over $L_2(X)$ and $L_2(Z)$, respectively. For $\theta\in\Theta = \{\theta | \sum_j \theta_j = 1, \theta_j\geq 0,\forall j\} $, we denote $h_{\theta} = \sum_j \theta_j h_j$. We use $e_j$ to denote the one-hot vector where that is zero except for the $j^{th}$ component, which equals to $1$. For a function class $\mathcal{F}$, we define the localized Rademacher complexity by $ \bar{R}_n (\delta; \mathcal{F}) := \mathbb{E}\big[ \mathbb{E}_\epsilon\big[\sup_{f\in \mathcal{F}, \|f\|_2 \leq \delta} \big|\frac{1}{n}\sum_{i=1}^n \epsilon_i f(x_i, z_i)\big|\big]\big],$ where $\epsilon_i$ are i.i.d. Rademacher random variables. For a function class $\mathcal{F}$ over $\mathcal{X}$ and $\mathcal{Z}$, we define its star hull by $\operatorname{star}(\mathcal{F}) = \{\gamma f, \gamma \in [0,1], f\in \mathcal{F} \}$. For a function class $\mathcal{F}$, we denote $\bar{\mathcal{F}} := \operatorname{star}(\mathcal{F} - \mathcal{F})$ to define its symmetrized star hull. We define the critical radius $\delta_{n,\mathcal{F}}$ of a function class $\mathcal{F}$ as any solution to the inequality $\delta^2 \geq \bar{R}_n(\operatorname{star}(\mathcal{F} - \mathcal{F}),\delta)$. We use $\mu$ to denote the Lebesgue measure.
As mentioned in \pref{sec:intro}, we aim to solve the following inverse problem with respect to $h$, known as the nonparametric IV regression:
While ${\mathcal{T}}$ and $r_0$ are unknown a priori, using i.i.d. observations $\{X_i, Y_i, Z_i\}_{i\in[n]}$, we aim to solve this equation. We denote its associated distributions by $g_0$, e.g., denote the conditional density of $X\in\mathcal{X}$ given $Z\in\mathcal{Z}$ by $g_0(x|z) \in \{\mathcal{X} \times \mathcal{Z} \rightarrow \mathbb{R}\}$. Throughout this work, we assume a solution to Equation (ref) exists.
Crucially, even though a solution to (ref) exists, it might not be unique. Hence, we propose to target a specific solution that achieves the least norm, defined as:
Note this least norm solution is well-defined, as it is defined by the projection of the origin onto a closed affine space $\mathcal{N}_{r_0}({\mathcal{T}}) \subset L_2(X)$. Indeed, with Assumption (ref), it is easy to prove that $h_0$ in (ref) always exists bennett2023minimax.
As we emphasize the challenges in \pref{sec:intro}, although there have been a lot of method that use minimax optimization for estimating $h_0$, when using general function approximation such as neural networks, the minimax optimization tends to be computationally hard lin2020near,jin2020local,lin2020gradient,diakonikolas2021efficient,razaviyayn2020nonconvex. Moreover, it remains unclear how to perform model selection for those methods. Hence, in this paper, we aim to propose a new method that can incorporate any function approximation for estimating the least square norm solution $h_0$ in (ref) with a strong convergence guarantee in $L_2(X)$ under mild assumptions (i.e., such as without the uniqueness of $h_0$) while allowing for model selection.
In this section, we introduce a two-stage algorithm, Regularized DeepIV (RDIV), aimed at obtaining the least square solution $h_0$ as defined in Equation (ref). Even though we borrow the DeepIV terminology from the prior work hartford2017deep, our method can be used with arbitrary function approximators and not necessarily neural network function spaces. Being inspired by the original constrained optimization (ref), we aim to solve a regularized version of the problem:
where $\mathcal{H} \subset L_2(X)$ represents a hypothesis class that consists of possible candidates for $h_0$, and $\alpha\in \mathbb{R}^+$ denotes a parameter controlling the strength of regularization. While this formulation itself has been known in the literature on general inverse problems cavalier2011inverse,mendelson2010regularization, we consider common scenarios in IV where both the conditional expectation operator ${\mathcal{T}}$ and the population expectation in Equation (ref) are unknown, and need to leverage dataset $\{X_i,Y_i,Z_i\}$.
To address this challenge, by integrating general function approximation such as neural networks, we propose a two-stage method, the Regularized Deep Instrumental Variable (RDIV), which is summarized in \pref{alg:dbiv-noniterative}. In the first stage, given a function class $\mathcal{G}$ comprising functions of the form $\big\{g: \mathcal{X}\times\mathcal{Z} \rightarrow \mathbb{R},\int_\mathcal{X} g(x|z) \mu(dx) =1 \text{ for all } z\big\}$, we aim to learn the conditional expectation operator ${\mathcal{T}}$ by estimating the ground-truth conditional density $g_0(x|z)$ from the dataset $\{X_i, Z_i\}_{i\in[n]}$ with MLE in Equation (ref). In the second stage, with the learned conditional density $\hat{g}$ in the first step, we learn $h_0$ by replacing expectation and ${\mathcal{T}}$ in Equation (ref) with empirical approximation and $\hat {\mathcal{T}}$, respectively, as shown in Equation (ref).
Importantly, our method does not necessitate a demanding computational oracle such as non-convex non-concave minimax or bilevel optimization, unlike many existing works for nonparametric IV with general function approximation lewis2018adversarial,xu2021deep, bennett2023minimax. Even when using neural networks for $\mathcal{G}$ and $\mathcal{H}$, we just need standard ERM oracles for density estimation or regression whose optimization is empirically known to be successful and theoretically more supported du2019gradient,chen2018convergence,zaheer2018adaptive,barakat2021convergence,wu2019global,zhou2018convergence,ward2020adagrad. We leave the numerical comparison between our method and existing NPIV methods hartford2017deep, dikkala2020minimax,xu2021deep,singh2019kernel in Appendix (ref).
In this section, we demonstrate a convergence result of our estimator $\hat{h}$ in RDIV to $h_0$ and derive its $L_2$ error rate after introducing several assumptions.
We commence by introducing the $\beta$-source condition, a concept commonly used in the literature on inverse problems carrasco2007linear,ito2014inverse,engl1996regularization,bennett2023source,liao2021instrumental, which mathematically captures the well-posedness of the function $h_0$.
In the following, we present its interpretation. First, as special cases, when $\mathcal{X},\mathcal{Z}$ are finite (e.g., discrete random variables), it holds when $\beta=\infty$. However, in our cases of interests where $\mathcal{X},\mathcal{Z}$ are not finite, this assumption restricts the smoothness of $h_0$. Intuitively, when the parameter $\beta$ is large, the function $h_0$ exhibits greater smoothness, and the assumption gets stronger, in the sense that eigenfunctions of $h_0$ relative to an operator ${\mathcal{T}}$ have smaller eigenvalues as explained in bennett2023minimax.
Next, we introduce another standard assumption as follows. This requires that the function classes $\mathcal{H}$ and $\mathcal{G}$ are well-specified. We will later consider misspecified cases as in \pref{sec:misspecified}.
The final assumption is as follows. This is standard in analyzing the convergence of nonparametric MLE wainwright2019high. We will later discuss how to relax such an assumption in Remark (ref) and Appendix (ref).
Finally, we present our guarantee for Algorithm (ref).
\paragraph{Sketch of the proof.} We now sketch the proof of Theorem (ref). Recall that $h_*$ is the optimizer of the Tikhonov-regularized loss function (ref). We first introduce the following lemma, which characterizes the bias caused by the regularization.
Therefore, recalling $\|\hat h- h_0 \|_2^2 \leq 2(\|h_* - h_0\|_2^2 + \|\hat{h} - h_*\|_2^2)$, we only need to bound $\|\hat{h} - h_*\|_2^2$. Utilizing the strong convexity of (ref), we have the following lemma:
Lemma (ref) shows we can bound $\|\hat{h} - h_{*}\|_2^2$ by two terms (a1) and (a2). Here (a1) is a centered empirical process, and (a2) is the error when estimating ${\mathcal{T}}$ by $\hat{\mathcal{T}}$. Utilizing localized concentration inequality and the boundedness of function class $\mathcal{H}$, we can bound (a1) by $O(\delta_{n,\mathcal{H}}^2 + \delta_{n,\mathcal{H}}\|\hat{h} - h_*\|_2)$. To control (a2), we prove the following lemma:
Note that such a bound is nontrivial since the standard analysis for MLE only results in a convergence rate in terms of Hellinger distance between the $\hat{g}$ and $g$, and does not directly provide a bound on $L_1$ norm of ${\mathcal{T}} h - \hat {\mathcal{T}} h$. With Lemma (ref), we can now bound (a2) from above with the critical radius of $\mathcal{G}$, and the $L_2$ distance between $\hat h$ and $h_*$. Finally, combining current arguments, Lemma (ref) and (ref), we have
for certain constants $c'$. By organizing the above equation, we have
Combine Lemma (ref) and Equation (ref), we further have $$ \textstyle \|\hat h - h_0\|_2^2 = O\bigg(\frac{\delta_n^2}{\alpha^2} + \alpha^{\min\{\beta,2\}}\bigg), $$ select $\alpha = \delta_n^{\frac{2}{2 + \min\{\beta, 2\}}}$ and we conclude the proof. $\null\nobreak\hfill\ensuremath{\blacksquare}$
The critical radius $\delta_n$ measures the statistical complexity of function class $\mathcal{H}$ and $\mathcal{G}$. For example, for parametric class or Gaussian Kernel, $\delta_n =\Tilde{O}(n^{-1/2})$, while for first order Sobolev class, $\delta_n = \tilde{O}(n^{-1/3})$ wainwright2019high, bartlett2002localized. In those cases, when $\beta \geq 2$, the final rate in $L_2$ metric will be $\tilde{O}(n^{-1/2})$ in the former case and $\tilde{O}(n^{-1/3})$ in the latter case, respectively. We now give the interpretation of our result. The bound of $\|\hat{h} - h_0\|_2^2$ consists of two terms. Term (i) comes from a statistical error to estimate $h_{*}$ from $\mathcal{H}$ and $\mathcal{G}$ (i.e., $\|\hat h- h_{*}\|_2$). Here, we use the strong convexity owing to Tikhonov regularization as it enables us to convert the population risk error to an error in $L_2$ metric as in Lemma (ref). Then, we properly bounded the population risk from above by the empirical process term properly as in Lemma (ref). While this $\delta^2_n$ rate is known as the standard fast rate in nonparametric regression wainwright2019high, our result is still non-trivial because we need to handle a statistical error term properly when approximating ${\mathcal{T}}$ with $\hat {\mathcal{T}}$, which comes from the MLE error in the form of Hellinger distance.
The term (ii) comes from the bias $\|h_0- h_{*}\|_2$ incurred by adding a Tikhonov regularization. This analysis has been used in existing works (e.g., cavalier2011inverse). Due to $\min(\beta,2)$, while we cannot leverage a high smoothness $\beta$ especially when $\beta \geq 2$, we will see how to leverage $\beta$ in such a case by introducing an iterative estimator in Section (ref).
We also compare our work to existing state-of-the-art convergence rate $O\big(\delta_n^{2\frac{\min\{\beta,1\}}{1+\min\{\beta,1\}}}\big)$ in bennett2023source, in which they employ a minimax-type algorithm. When $\beta \geq 2$, we achieve the same rate. We also remark that although our rate is slightly slower than theirs when $\beta \leq 2$, our method does not require a minimax-optimization oracle and can be incorporated with method selection methods. Besides, we will show that our method can achieve a state-of-the-art rate in our extension to iterative estimator in \pref{sec: iterative}.
Next, we establish the finite sample result when Assumption (ref) does not hold, i.e., function classes $\mathcal{H}$ and $\mathcal{G}$ are misspecified. This result serves as an important role in formalizing the model selection procedure in \pref{sec:model-selection}.
The bound for $\|\hat{h} - h_0\|^2_2$ consists of three terms: term (b1) measures the statistical deviation of a normalized empirical process, term (b2) measures the regularization error caused by Tikhonov regularization and term (b3) measures the effect of model misspecification. Here term (b3) has a poly$(\frac{1}{\alpha})$ dependency. This is because model misspecification causes a higher population risk in both stage 1 and 2 of Algorithm (ref). Hence, the more convex the loss function, the lesser the shift in the optimizer. The readers may notice that term (b2) is slightly slower than the original bias term in Theorem (ref). This is because the difference of the optimal value in (ref) due to misspecification of $\mathcal{H}$ is of order $O(\alpha^{\min\{\beta+1,2\}}+ \epsilon_{\mathcal{H}}^2) $, as we will show in Lemma (ref) in the Appendix. By the $\alpha$-strong convexity endowed by Tikhonov regularization, this results in a shift of $h_*$ of magnitude $O\big(\alpha^{\min\{\beta+1,2\}-1}+ \epsilon_{\mathcal{H}}^2/\alpha\big)$.
Theorem (ref) is particularly useful when we apply estimators based on sample-dependent function classes $\mathcal{H}$ and $\mathcal{G}$ (a.k.a. sieve estimators) that approximate certain function spaces. For example, $\mathcal{H}$ can be linear models with polynomial basis functions that take the form $\langle \phi(X), \theta\rangle$, which can gradually approach H\"{o}lder or Sobolev balls, and $\mathcal{G}$ can be a set of neural networks with a growing dimension chen2007large,chen2022nonparametric,schmidt2020nonparametric. More specifically, when $X$ and $Z$ are bounded, and $h_0$ and $g_0$ are $s$-H{\"o}lder smooth, it is well known that a deep ReLU neural network with depth $O(\log(1/\epsilon))$, width $O(d\epsilon^{-d/s})$ and weights bounded by $\tilde{O}(1)$ could satisfy the approximation error in Theorem (ref) schmidt2019deep, recall that $d$ is the dimension of $X$ and $Z$ . In that case, $\delta_n^2 = \tilde{O}(\epsilon^{-d/s}/n)$ bartlett2002localized, chen2022nonparametric. Choosing the architecture of the neural network according to $\epsilon = \tilde{O}(n^{-1/(1+d/s)})$, then Theorem (ref) shows that by setting $\alpha = O(n^{\frac{1}{(1+d/\alpha)(\min\{\beta+1,2\}+1)}})$, we have $\|\hat{h} - h_0\|_2^2 = \tilde{O}(n^{\frac{\min\{\beta+1,2\}-1}{(1+d/s)(\min\{\beta+1,2\}+1)}})$.
One advantage of employing the proposed two-staged algorithm is that it enables model selection, which is not attainable when a minimax approach is used. In this section, we explain how we perform model selection. We focus on the model selection for the second stage, as the conditional density $\hat{g}$ from the first stage can be selected via existing methods for model selection for maximum likelihood estimators (e.g. MLEselection,conditiondensitySelection,vijaykumar2021localization).
With an MLE-based estimator $\hat{g}$ obtained from the first stage in Algorithm (ref), we consider model selection using the regularized loss in the second stage, with theoretical guarantees in the $\|\cdot\|_{2}$ metric. More concretely, given a choice of $M$ candidate models $\{h_1, \dots, h_M\}$ and a validation dataset $\{X'_i,Y'_i,Z'_i\}_{i=1}^n$ (distinct from the one used for training models $\{h_i\}$ and $\hat{g}$), the goal is for the final output of the model selection algorithm to achieve oracle rates with respect to the minimal misspecification error.
We present our algorithm in Algorithm (ref). We provide two options for model selection: Best-ERM and Convex-ERM. Best-ERM selects the model that minimizes the regularized loss on a validation set, while Convex-ERM constructs a convex aggregate of the candidate models that minimizes the regularized loss on a validation set.
We explain its implications. Most importantly, our obtained rate is the best (i.e., oracle rate) among rates when invoking a result of (convergence result for RDIV in \pref{thm: mle-misspec} with misspecified model) for each function class $\mathcal{H}_i$. Some astute readers might wonder whether we can just invoke \pref{thm: mle-misspec} by making new function classes $\mathcal{H}_{\text{best}}:=\{h_{\theta}:\theta=e_1,\dots,e_M\}$ or $\mathcal{H}_{\text{conv}}:=\{h_{\theta}:\sum_j \theta_j = 1, \theta_j \geq 0 \}$, and bound the misspecification error $\epsilon_{\mathcal{H}_{\text{conv}}}$ or $\epsilon_{\mathcal{H}_{\text{best}}}$ by $\|h_j - h_0\|$ will lead to a slower rate with an extra factor of $\frac{1}{\alpha}$. The key is only to handle the misspecification error once to avoid the $\frac{1}{\alpha}$ factor by deferring the invocation of strong convexity and working with the excess risk (difference in the expected loss) instead of the $L_2$ difference.
One drawback of the result so far is its lack of adaptability to the degree of ill-posedness in the inverse problem, especially for larger values of $\beta$ corresponding to milder problems, when $\beta \geq 2$. To address this issue, in this section, we further generalize our results in Section (ref) and (ref), and propose an iterated Regularized Deep method, which is summarized in Algorithm (ref). In this algorithm, instead of targeting (ref), we target $h_{m, *}$, which is given by the following recursive least square regression with Tikhonov regularization:
and we set $h_{-1,*} = 0$. This is the recursive version of the previous regularized objective in Equation (ref), by using Tikhonov regularization around a prior target $h_{m-1,*}$ instead of $0$. Then, with the learned conditional density $\hat{g}$ by MLE in Equation (ref), we construct an estimator in (ref) by replacing expectation and an operator ${\mathcal{T}}$ with empirical approximation and the learned operator $\hat {\mathcal{T}}$, respectively, in Equation (ref).
Now, we delve into estimating the finite sample convergence rate of Algorithm (ref). Our findings are summarized in the following theorem.
Importantly, we can have a rate $O\big(\delta_n^{\frac{2\beta}{2 + \beta}}\big)$ in relatively mild conditions while the previous Theorem (ref) (non-iteratie version) can only allow for $O\big(\delta_n^{\frac{2\min(\beta,2)}{2 + 2\min(\beta,2)}}\big)$, and cannot fully leverage the well-posedness of $h_0$, illustrated by the source condtion $\beta$. Indeed, if we choose the iteration number $m = \lceil \min\{\beta/ 2, \log\log(1/\delta_n)\} \rceil$, then we get a rate of $$\|\hat{h}_m - h_0\|_2^2 = O\bigg(\min\{16^\beta, \log (1/\delta_n)\}\delta_n^{\frac{2\min\{\beta, 2m\}}{2 + \min\{\beta, 2m\}}}\bigg) .$$ Hence for any constant $\beta$, as $n$ grows, eventually $\log\log 1/\delta_n \geq \beta $, and we get the rate of $O\big(\delta_n^{\frac{2\beta}{2 + \beta}}\big)$. This rate can be achieved even if $\beta$ grows with $n$, as long as it grows slower than $O(\log\log 1/\delta_n)$. If $\delta_n = O(n^{-\iota})$ for some $\iota > 0 $, e.g. RKHS or first order Sobolev space wainwright2019high, then we note that we can set $m = \lceil \min\{\beta/ 2, \sqrt{\log(1/\delta_n)}\} \rceil$, and $16^{\sqrt{\log(1/\delta_n)}} = O(n^{\epsilon})$ for any $\epsilon > 0$, thus we still obtain a rate of $O(\delta_n^{\frac{2\beta}{2+\beta}})$ when $\sqrt{\log(1/\delta_n)} \geq \beta/2$. In such a case, we can obtain a $O(\delta_n^{\frac{2\beta}{2+\beta}})$ rate even $\beta$ grows with $n$, as long as it grows slower than $\sqrt{\log(1/\delta_n)}$.
Our results for the iterative estimator match the state-of-the-art convergence rate with respect to $L_2$ norm for an iterative estimator in bennett2023source. However, their method requires a minimax computation oracle, while our method does not.
In this section, we evaluate our proposal by numerical simulation. In particular, we present the performance of RDIV when we use neural networks as the function approximator and the validity of the proposed model selection procedure. We show that with model selection, our method can achieve state-of-the-art performance in a wide range of data-generating processes.
\paragraph{Experiment Design.} In our experiment, we test our method on a synthetic dataset. We adjust the data generating process (DGP) for proximal causal inference used in cui2020semiparametric,miao2018confounding,deaner2021many. Concretely, we generate multi-dimensional variables $U',S', W',Q',A$, where $U$ is an unobserved confounder, $S'\in d_{S}$ is the observed covariate, $W'\in d_{W}$ is the negative control outcomes, $Q'\in d_{Q}$ is the negative control actions, and $A$ is the selected treatment, as described in \pref{fig:causal_dag}. We left the detailed generation process in Appendix (ref). For a detailed understanding of this setup, we refer the reader to Section 2 of kallus2021causal. It is well known that there exists a bridge function $h_0'$ such that the following moment condition holds cui2020semiparametric,kallus2021causal: $$ \mathbb{E}[Y - h_0'(W',A,S')|Q',A,S'] = 0, $$ which allows the concrete form of (ref). To introduce nonlinearity, we transform $(S', W',Q')$ into $(S,W,Q)$ via $S = g(S'), W = g(W'), Q=g(Q')$, where $g(\cdot)$ is a nonlinear invertible function applied elementwise to $S',W',Q'$ respectively. We consider several forms of $g(\cdot)$, including identity, polynomial, sigmoid design, and exponential function. In the final data, we only observe $(S,W,Q)$ but not $(S',W',Q')$. Here we use 6 different $g(\cdot)$: Id$(t) = t$, Poly$(t)$ = $t^3$, LogSigmoid$(t)= \log(1+ |16*x - 8|) \cdot \operatorname{sign}(x)$, Piecewise$(t) = 3(x-2)1_{x\leq1}+ \log(8x-8)1_{x\geq1}$, Sigmoid$(t) = \frac{5}{1+\exp(-0.1*x)}$ and CubicRoot = $x^{1/3}$.
\paragraph{Methods to compare.}
In this experiment, our goal is to estimate the counterfactual mean parameter $\mathbb{E}[Y(1)]$, which is unique as long as (ref) holds. We learn $h_0$ in (ref) by RDIV, which corresponds to the procedure in Algorithm (ref) with MLE for conditional density estimation. We show results for different values for $\alpha \in \{0.01,0.1\}$, and compare the performance of our approach to that of several different methods, including KernelIV singh2019kernel, DeepIV hartford2017deep, DeepFeatureIV xu2021deep, and AGMM DikkalaNishanth2020MEoC. Note that DeepIV can be viewed as a special case of our methods, with $\alpha$ fixed to be 0. In the first stage of our algorithm, we use a three-layer mixture density network hartford2017deep,rothfuss2019conditional as the approximator of the conditional density. In the second stage, we use a three-layer fully-connected neural network as the approximators for RDIV, DeepIV, AGMM, and DFIV. We present the results of our method and its comparison with previous benchmarks in terms of MSE normalized by the true estimand value in Table (ref)-(ref). Every estimate is calculated by 100 random replications. The confidence interval is calculated by 2 times the standard deviation.
\paragraph{Hyperparameter settings.} For RDIV, we use Adam as the optimizer for both density estimation and Tikhonov regression, with a default learning rate of $10^{-4}$, a batch size of $50$, and a training epoch of $300$. We will show how to choose these hyperparameters with our model selection procedure (Algorithm (ref)) in Section (ref). For all baselines except for AGMM, we adapt the hyperparameters in their original codebase. For AGMM, we tune the learning rate for the learner and adversary for every $g(\cdot)$ independently. We follow singh2019kernel to use Gaussian RKHS for function approximation and their method for tuning the regularization parameter. When $n=500$, the learning rate of the learner and adversary in AGMM are manually set to $10^{-4}$ for LogSigmoid, Piecewise, and Sigmoid, and $10^{-3}$ for Id, Poly, and CubicRoot. When $n=1000$, the learning rate of the learner and adversary in AGMM are manually set to $10^{-4}$ for Piecewise and Sigmoid, and $10^{-3}$ for LogSigmoid, Piecewise, and CubicRoot. The training parameter of DFIV is adopted from xu2021deep. Note that tuning DFIV is highly intractable in practice, as their method is essentially a bilevel optimization, which is known to be hard to solve hong2023two.
First, we can observe that although our estimator resembles DeepIV, the later fix $\alpha =0$ in (ref), RDIV outperforms DeepIV for all $g(\cdot)$. This is due to the nonzero regularization term, which improves the performance of our estimator by a better tradeoff between bias and variance. Second, in most cases, AGMM and DFIV are outperformed by algorithms that only need single-level optimization (RDIV, KernelIV, DeepIV). This would be because, in these methods, optimization of the loss function is much harder, which results in the inaccuracy of estimators. Thirdly, while it is seen that Kernel IV is comparable to RDIV in some scenarios such as in Table (ref) and (ref), in the next section, we will show that our RDIV equipped with a model selection procedure can generally outperform KernelIV.
We also report our results in model selection for the second stage by implementing Best-ERM in Algorithm (ref) and demonstrate how it improves our results. Specifically, our models $h_1, \dots, h_M$ are trained by different hyperparameters. First, we employ model selection for the density function by Best ERM. Then with the trained density function in the first stage, we further apply Best ERM to the models in the second stage. In the model selection experiments, we fix the dimension of our dataset to be $d_S = d_Q = 20$, $d_W = 10$. We compute the mean and confidence interval with 10 independent trials. We set the candidate training parameters as follows: the number of epochs $\in\{300,400\}$, the batch size for the 1st stage $\in\{30,50\}$ and the batch size for the 2nd stage $\in\{50,60,100\}$, the learning rate $\in\{10^{-4},10^{-3}\}$, the number of mixture components $\in\{40,50,60\}$. As shown in Table (ref), when RDIV is equipped with model selection techniques, our method outperforms KernelIV in all but one case when the dataset size is 500, and outperforms KernelIV in 3 out of 6 settings when the dataset size is 1000. Our approach demonstrates its effectiveness by outperforming previous benchmarks across a diverse set of Data Generating Processes (DGP). This achievement is attributed to both the ease of optimization of RDIV and its theoretically sound integration with model selection procedures.