EconBase
← Back to paper

Stochastic Optimization Algorithms for Instrumental Variable Regression with Streaming Data

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.

51,767 characters · 6 sections · 71 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Stochastic Optimization Algorithms for Instrumental Variable Regression with Streaming Data

\footnotetext{XC and AR contributed equally to this work.}

abstractWe develop and analyze algorithms for instrumental variable regression by viewing the problem as a conditional stochastic optimization problem. In the context of least-squares instrumental variable regression, our algorithms neither require matrix inversions nor mini-batches and provides a fully online approach for performing instrumental variable regression with streaming data. When the true model is linear, we derive rates of convergence in expectation, that are of order $\mathcal{O}(\log T/T)$ and $\mathcal{O}(1/T^{1-\iota})$ for any $\iota>0$, respectively under the availability of two-sample and one-sample oracles, respectively, where $T$ is the number of iterations. Importantly, under the availability of the two-sample oracle, our procedure avoids explicitly modeling and estimating the relationship between confounder and the instrumental variables, demonstrating the benefit of the proposed approach over recent works based on reformulating the problem as minimax optimization problems. Numerical experiments are provided to corroborate the theoretical results.

Introduction

Instrumental variable analysis is widely used in fields like econometrics, health care, social science, and online advertisement to estimate the causal effect of a random variable, $X$, on an outcome variable, $Y$, when an unobservable confounder influences both. By identifying an instrumental variable correlated with the variable $X$ but unrelated to the confounders, researchers can isolate the exogenous variation in $X$ and estimate a causal relationship between $X$ and $Y$. In the context of regression, Instrumental Variable Regression (IVaR) addresses endogeneity issues when an independent variable is correlated with the error term in the regression model, leveraging an instrument variable $Z$ such that $Y$ is independent of $X|Z$. In this paper, we focus on the following statistical model:

align[align omitted — 113 chars of source]

where $X \in \mathbb{R}^{d_x}$ and $\epsilon_1$ are correlated and $\epsilon_2$ is a centered unobserved noise (independent of $Z\in\mathbb{R}^{d_z}$), leading to confounding in the model between $X$ and $Y \in \mathbb{R}$. Here $\epsilon_1$ and $\epsilon_2$ are dependent, and $\theta^*$ and $\gamma^*$ are true parameters for the respective function $g$ and $h$. Our goal is to design efficient algorithms that recovers $\theta^*$ from the data.

Traditionally, IVaR algorithms are based on two-stage estimation procedures, where we first regress $Z$ and $X$ to obtain an estimator $\hat X$, and then regress $\hat X$ and $Y$, with the essence that $\hat X$ is independent of $Y$, and thus eliminating the aforementioned endogeneity of the unknown confounder. A vast literature has devoted to understanding the two-stage approaches hall2005nonparametric, darolles2011nonparametric, hartford2017deep, with the parametric two-stage least-squares (2SLS) procedure being the most canonical one angrist1995two. The main drawback of this approach is that the second-stage regression problem is affected by the estimation error from the regression problem corresponding to first stage. In fact, angrist2009mostly call the first stage regression as “forbidden regression”, due to the concerns in estimating a nuisance parameter.

Considering the squared loss function, muandet2020dual formulate the IVaR problem as a conditional stochastic optimization problem hu2020biased:

align[align omitted — 125 chars of source]

However, muandet2020dual did not solve problem (ref) efficiently, and resort to reformulating (ref) further as a minimax optimization problem. Indeed, they mention explicitly in their work that “it remains cumbersome to solve (ref) directly because of the inner expectation”. Then, they leverage the Fenchel conjugate of the squared loss,

leading to a minimax optimization with maximization over a continuous functional space. Following dai2017learning, muandet2020dual propose to use reproducing kernel Hilbert space (RKHS) to handle the maximization over continuous functional space. See also lewis2018adversarial,bennett2019deep,dikkala2020minimax,liao2020provably,bennett2023minimax for similar minimax approaches. The issue with such an approach is that approximating the dual variable via maximization over continuous functional space inevitably introduces approximation error. Hence, although there is no explicit nuisance parameter estimation step like in the two-stage approach, there is an implicit one, which makes the minimax approach less appealing as an alternate to the two-stage procedures.

In this work, contrary to the claim made in muandet2020dual that problem (ref) is cumbersome to solve, we design and analyze efficient streaming algorithms to directly solve the conditional stochastic optimization problem in (ref). Direct application of methods from hu2020biased for solving (ref) is possible, yet their approach utilizes nested sampling, i.e., for each sample of $Z$, hu2020biased generate a batch of samples of $X$ from $\PP(Z|X)$, to reduce the bias in estimating the composition of non-linear loss function with conditional expectations. Thus their methods are not suitable for the streaming setting that we are interested in. Considering (ref), we first parameterize the function class $\mathcal{G}:=\{g(\theta;X) \mid \theta \in\RR^{d_\theta}\}$. Now, defining $F(g)\coloneqq F(\theta)$, we observe that the gradient $\nabla F(\theta)$ admits the following form

align[align omitted — 154 chars of source]

which implies that one does not need the nested sampling technique to reduce the bias. However, the presence of product of two conditional expectations $\EE_{X\mid Z}[g(\theta;X)]$ still causes significant challenges in developing stochastic estimators of the above gradient in the streaming setting. In this work, we overcome this challenge and develop two algorithms that are applicable to the streaming data setting avoiding the need for generating batches of samples of $X$ from $\PP(Z|X)$.\\

\noindentContributions. We make the following contributions in this work.

itemize[leftmargin=0.15in] • Two-sample oracles: Our first algorithm leverages the observation that if we have access to a two-sample oracle that outputs two samples $X$ and $X'$ that are independent conditioned on the instrument $Z$, we can immediately construct an unbiased stochastic gradient estimator of the gradient in (ref). Based on this crucial observation, we propose the Two-Sample One-stage Stochastic Gradient IVaR (TOSG-IVaR) method (Algorithm (ref)) that avoids explicitly having to estimate or model the relationship between $Z$ and $X$ thereby overcoming the “forbidden regression” problem.. Under standard statistical model assumptions, for the case when $g$ is a linear model, we establish rates of convergence of order $\mathcal{O}(\log T/T)$ for the proposed method, where $T$ is the overall number of iterations; see Theorem (ref). • One-sample oracles: In the case when we do not have the aforementioned two-sample oracle, we estimate the stochastic gradient in (ref) by using the streaming data to estimate one of the conditional expectations, and the corresponding prediction to estimate the other, resulting in the One-Sample Two-stage Stochastic Gradient IVaR (OTSG-IVaR) method (Algorithm (ref)). Assuming further that the $X$ depends linearly on the instrument $Z$, we establish a rate of convergence of order $\mathcal{O}(1/T^{1-\iota})$, for any $\iota>0$; see Theorem (ref).

Literature Review

\noindentIVaR analysis. Instrumental variable analysis has a long history, starting from the early works by wright1928tariff and reiersol1945confluence. Several works considered the aforementioned two-stage procedure for IVaR; a summary could be found in the work by angrist2009mostly. Nonparametric approaches based on wavelets, splines, reproducing kernels and deep neural networks could be found, for example, in the works by hartford2017deep,singh2019kernel,bennett2019deep, muandet2020dual,mastouri2021proximal,xu2021learning,zhu2022causal,peixoto2024nonparametric. Another popular approach for IVaR is via Generalized Method of Moments (GMM); see, for example, chen2012estimation,bennett2019deep,dikkala2020minimax for an overview. Such approaches essentially reformulate the problem as a minimax problem and hence suffer from the aforementioned “forbidden regression” problem.\% similar to muandet2020dual, i.e., the “forbidden regression” problem.

Identifiability conditions for IVaR. Several works in the literature have also focused on establishing the identifiability conditions for IVaR in the parametric and the nonparametric setting. Regardless of the procedure used, they are invariably based on certain source conditions motivated by the inverse problems literature (see, for example, carrasco2007linear,chen2011rate, bennett2023minimax) or the related problem of completeness conditions, which posits that the conditional expectation operator is one-to-one babii2017completeness,liao2020provably. Semi-parametric identifiability is also considered recently in the work of cui2023semiparametric. Our focus in this work is not focused on the identifiability; for the formulation (ref) that we consider, muandet2020dual provide necessary conditions for identifiability that we adopt.\\

\noindentStochastic optimization with nested expectations. Recently, much attention in the stochastic optimization literature has focused on optimizing a nested composition of $T$ expectation functions. Sample average approximation algorithms in this context are considered in the works of ermoliev2013sample and hu2020sample. Optimal iterative stochastic optimization algorithms for the case of $T=2$ were by derived by ghadimi2020single. For the general $T\geq 1$ case, wang2017stochastic provided sub-optimal rates, whereas balasubramanian2022stochastic derived optimal rates; see also zhang2021multilevel and chen2021solving for related works under stronger assumptions, and ruszczynski2021stochastic for similar asymptotic results. While the above works required certain independence assumptions regarding the randomness across the different compositions, hu2020biased,hu2024contextual studied the case of $T=2$ where the the randomness are generically dependent. They termed this problem setting as conditional stochastic optimization, which is the framework that the IVaR problem in (ref) falls in. Compared to prior works, for e.g., ghadimi2020single and balasubramanian2022stochastic, in order to handle the dependency between the levels, hu2020biased require mini-batches in each iteration, making their algorithm not immediately applicable to the purely streaming setting. In this work, we show that despite the problem (ref) being a conditional stochastic optimization problem, mini-batches are not required due the additional favorable quadratic structure available in IVaR.\\

Streaming IVaR. venkatraman2016online,della2023online analyzed streaming versions of 2SLS in the online\footnote{Their notion of online is from the literature on online learning shalev2012online.} and adversarial settings. Focusing on linear models, venkatraman2016online provide preliminary asymptotic analysis assuming access to efficient no-regret learners, while della2023online provide regret bounds under the strong assumption that the instrument is almost surely bounded. Furthermore, our algorithms have significantly improved per-iteration and memory complexity compared to della2023online; see Sections (ref) and (ref) for details. chen2023sgmm developed stochastic optimization algorithms for the GMM formulation and provide asymptotic analysis. Their algorithm requires access to an offline dataset for initialization and is hence not fully online. The above works (i) do not focus on avoiding the forbidden regression problem and (ii) do not view IVaR via the conditional stochastic optimization lens, like we do.

Two-sample One-stage Stochastic Gradient Method for IVaR

Recall that our goal is to solve the objective function given in (ref). By muandet2020dual, the optimal solution of (ref) gives the true underlying causal relationship under the following assumption.

assumption(Identifiability Assumption) \begin{itemize}[noitemsep,leftmargin=0.12in] • The conditional distribution $\PP_{Z\mid X}$ is continuous in $Z$ for any value of $X$. • The function class $\mathcal{G}:=\{g(\theta;X) \mid \theta \in\RR^{d_\theta}\}$ is correctly specified, i.e., it includes the true underlying relationship between $X$ and $Y$. \end{itemize}

Notice that both assumptions are standard in the IVaR literature newey2003instrumental,chen2012estimation,muandet2020dual, and makes the objective in (ref) is the meaningful for IVaR. However, muandet2020dual resort to reformulating the objective function in (ref) as a minimax optimization problem as described in Section (ref). While their original motivation was to avoid two-state estimation procedure and avoid the “forbidden regression”, their minimax reformulation ends up having to solve a complicated approximation of the original objective resulting in having to characterize the approximation error which is non-trivial.\\

Algorithm and Analysis. Our aim in this work is to directly solve the original problem in (ref), leveraging the structure provided by the quadratic loss. Given the gradient formulation in (ref), a natural way to build unbiased gradient estimator is to generate $X$ and $X^\prime$, two independent samples of $X$ from the conditional distributions $\PP_{X\mid Z}$, for a given realization of $Z$ and generate one sample of $Y$ from the conditional distribution $\PP_{Y\mid X}$. Then, an unbiased gradient estimator is

equation[equation omitted — 108 chars of source]

This could be plugged into the standard stochasic gradient descent algorithm, which give us the Two-sample Stochastic Gradient Method for IVR (TSG-IVaR) method illustrated in Algorithm (ref). In particular, the algorithm never requires estimating (or modeling) the relationship between $X$ and $Z$ as needed in the two-stage procedure angrist2009mostly and the minimax formulation based procedures muandet2020dual,lewis2018adversarial,bennett2019deep,dikkala2020minimax,liao2020provably,bennett2023minimax. Furthermore, this viewpoint not only provides a novel algorithm for performing IV regression, but also provides a novel data collection mechanism for the practical implementation of IVaR. In addition, such a two-sample gradient method is not very restrictive when the instrumental variable $Z$ takes value in a discrete set. In this case, to implement the two-sample oracle, it is enough simply pick two sets of samples $(X,Y,Z)$ and $(X',Y',Z)$ for which $Z$ has repeated observations (which is possible when $Z$ is a discrete random variable) from a pre-collected dataset. To demonstrate the convergence rate of Algorithm (ref), we first consider the case when $g$ is a linear function, i.e., $g(\theta;X)=X^\top \theta$. We make the following assumptions.

algorithm[algorithm omitted — 672 chars of source]
assumptionSuppose there exists $\mu>0$ such that $ \EE_Z\Big[\EE_{X\mid Z}[X]\cdot \EE_{X\mid Z}[X]^\top\Big] \succeq \mu I.$
assumptionLet $(\vartheta_1, \vartheta_2, \vartheta_3, \vartheta_4) \in \mathbb{R}_+^4$. For any $Z$, $X'$ and $X$ i.i.d. generated from $\PP_{Z\mid X}$ , and $Y$ generated from $\PP_{Y\mid X}$, and for constants $C_{x},C_y,C_{xx},C_{yx}>0$, we have \begin{align} &\EE\Big[\norm{X' X^\top - \EE_{X|Z}[X]\EE_{X|Z}[X]^\top}^2\Big]\leq C_{x}d_x^{\vartheta_1}, \\ &\EE\Big[\norm{YX - \EE_{Y|Z}[Y]\EE_{X|Z}[X]}^2\Big]\leq C_y d_x^{\vartheta_2}, \\ &\EE\Big[\norm{\EE_{X\mid Z}[X]\cdot \EE_{X\mid Z}[X]^\top -\EE_Z\Big[\EE_{X\mid Z}[X]\cdot \EE_{X\mid Z}[X]^\top\Big]}^2\Big] \leq C_{xx}d_z^{\vartheta_3}, \\ &\EE\Big[\norm{\EE_{Y\mid Z}[Y]\cdot \EE_{X\mid Z}[X] -\EE_Z\Big[\EE_{Y\mid Z}[Y]\cdot \EE_{X\mid Z}[X]\Big]}^2\Big] \leq C_{yx}d_z^{\vartheta_4}, \end{align} where $\|\cdot\|$ denotes the Euclidean norm and operator norm for a vector and matrix respectively.

The above assumptions are mild moment assumptions required on the involved random variables. The following result demonstrates that Assumptions (ref) and (ref) are naturally satisfied under even under non-linear modeling assumption on (ref). We defer its proof to Section (ref).

lemmaSuppose there exist $\theta_*\in \RR^{d_x},\ \gamma_*\in\RR^{d_z\times d_x},$ a non-linear map $\phi: \RR^{d_x}\rightarrow \RR^{d_x}$, and a positive semi-definite matrix $\Sigma\in \RR^{d_z\times d_z}$ such that \begin{align} &\EE_Z\Big[\phi(\gamma_*^\top Z)\cdot \phi(\gamma_*^\top Z)^\top\Big]\succeq \mu I,\ \EE[\norm{\phi(\gamma_*^\top Z)}^2] = \cO(d_x), \notag \\ &Z \sim \cN(0, \Sigma),\ X = \phi(\gamma_*^\top Z) + \epsilon_2,\ Y = \theta_*^\top X + \epsilon_1,\ \epsilon_2\sim \cN(0, \sigma_{\epsilon_2}^2I_{d_x}),\ \epsilon_1 \sim \cN(0, \sigma_{\epsilon_1}^2), \end{align} where $\epsilon_1, \epsilon_2$ are independent of $Z$ and \begin{equation} \EE\left[\epsilon_1^2\norm{\epsilon_2}^2\right]\leq \sigma_{\epsilon_1,\epsilon_2}^2d_x,\ \EE\left[\norm{\phi(\gamma_*^\top Z)\cdot \phi(\gamma_*^\top Z)^\top - \EE[\phi(\gamma_*^\top Z)\cdot \phi(\gamma_*^\top Z)^\top]}^2\right]\leq Cd_z, \end{equation} then Assumptions (ref) and (ref) hold with $\vartheta_1=\vartheta_2=2$ and $\vartheta_3=\vartheta_4=1$. if $\phi$ is an identity map, then the conditions involving $\phi$ become $ \gamma_*^\top \Sigma\gamma_*\succeq \mu I,\ \tr(\gamma_*^\top\Sigma\gamma_*) = \cO(d_x),\ \EE\left[\norm{ZZ^\top - \Sigma}^2\right]\leq Cd_z.$
assumptionThe tuple $(Z_t, X_t, X'_t, Y_t)$ is independent and identically distributed, across $t$.

The above assumption is standard in the stochastic approximation, statistics and econometrics literature. It could be further relaxed to Markovian-type dependency assumptions, following techniques in the works of duchi2012ergodic,sun2018markov,even2023stochastic, roy2022constrained; we leave a detailed examination of the Markovian streaming setup as future work. Under the above assumptions, we have the following result demonstrating the last-iterate global convergence of Algorithm (ref).

theoremSuppose Assumptions (ref), (ref), and (ref) hold. In Algorithm (ref), defining $\sigma_1^2 \coloneqq 2C_{x}d_x^{\vartheta_1} + 2C_{xx}d_z^{\vartheta_3}$ and $\sigma_2^2 \coloneqq C_y d_x^{\vartheta_2} + C_{yx}d_z^{\vartheta_4}$, set $\alpha_t \equiv \alpha = \frac{\log T}{\mu T}\leq \frac{\mu}{\mu^2 + 3\sigma_1^2}$. Then, we have \begin{align*} \EE\big[\norm{\theta_T - \theta_*}^2\big]\leq \frac{\EE\big[\norm{\theta_0 - \theta_*}^2\big]}{T} + \frac{3\norm{\theta_*}^2(\sigma_1^2 + \sigma_2^2)\log T}{\mu^2T}. \end{align*}

\noindentProof techniques. In the analysis of Theorem (ref), the following decomposition (see (ref) for the derivation) plays a crucial role:

align*[align* omitted — 463 chars of source]

where $A_t$ corresponds to deterministic component, and $B_t$ corresponds to the stochastic component arising due to the use of stochastic gradients. Standard assumptions on the variance of the stochastic gradient made in the stochastic optimization literature include the uniformly bounded variance assumption lan2020first and the expected smoothness condition khaled2020better. In the IVaR setup, such standard assumptions do not hold as $\theta_t$ potentially can be unbounded and thus the gradient estimator can be unbounded. Hence, we establish our results under natural statistical assumptions arising in the context of the IVaR problem, which form the main novelty in our analysis. Furthermore, compared to muandet2020dual, notice that we use two samples of $X$ from the conditional distribution $\PP_{X|Z}$ and achieve an $\tilde \cO(1/T)$ last iterate convergence rate to the global optimal solution, which is the true underlying causal relationship under Assumption (ref). In comparison, muandet2020dual only provide asymptotic convergence result to the optimal solution of an approximation problem.\\

Additional discussion. It is interesting to explore other losses beyond squared loss (for example to handle classification setting centorrino2021nonparametric), potentially using the Multilevel Monte Carlo (MLMC) based stochastic gradient estimators. While hu2021bias, develops such algorithms, the main challenge is about how to avoid mini-batches required in their work leveraging the problem structure in instrumental variable analysis. Furthermore, in the case when $g(\theta;X)$ is parametrized by a non-linear models, for instance, a neural network, we provide local convergence guarantees under additional stronger conditions made typically in the stochastic optimization literature.

assumptionLet the following assumptions hold: \begin{itemize}[noitemsep, leftmargin=0.21in] • Function $F(\theta)$ is $\ell$-smooth. • The iterates $\{\theta_t\}_{t=1}^{T+1}$ generated by Algorithm (ref) are in a compact set $A$. • The random objects $X|Z$ and $Y|Z$ have bounded variance for any $Z$, i.e., there exist $\sigma > 0$ such that \begin{align*} \EE\left[\norm{X - \EE\left[X\mid Z\right]}^2\mid Z\right]\leq \sigma^2,\ \EE\left[\norm{Y - \EE\left[Y\mid Z\right]}^2\mid Z\right]\leq \sigma^2. \end{align*} \end{itemize}
propositionSuppose Assumptions (ref), (ref), and (ref) hold. Choosing $\alpha_t\equiv \alpha = \mathcal{O}\left(\frac{1}{\sqrt{T}}\right)$, for Algorithm (ref) we have \begin{align*} \min_{1\leq t\leq T}\EE\left[\norm{\nabla F(\theta_t)}^2\right] = \mathcal{O}\left(\frac{1}{\sqrt{T}}\right). \end{align*}

The proof of the proposition is immediate. Note that under Assumption (ref), we can deduce that the unbiased gradient estimator $v(\theta) = (g(\theta; X) - Y)\nabla_{\theta}g(\theta; X')$ has a bounded variance since

align*[align* omitted — 316 chars of source]

where the variance and expectation are taken conditioning on $Z$ and $\theta$, and $\sigma_v>0$ is a constant that only depends on $\sigma$, function $g$ and the compact set $A$ in Assumption (ref). Then one can directly follow the analysis of non-convex stochastic optimization (see, for example, ghadimi2013stochastic) to obtain Proposition (ref). Relaxing the Assumption (ref) (typically made in the stochastic optimization literature) with more natural assumptions on the statistical model and obtaining a result as in Theorem (ref) for the non-convex setting is left as future work.

One-sample Two-stage Stochastic Gradient Method for IVaR

We now examine designing streaming IVaR algorithm with access to the classical one-sample oracle, i.e., we observe a streaming set of samples $(X_t,Y_t,Z_t)$ at each time point $t$. Note that in this case, using the same $X_t$ (instead of $X'_t$) in (ref) makes the stochastic gradient estimator biased.\\

Intuition. Consider the case of linear models, i.e., $Y=\theta_*^\top X +\epsilon_1$ with $X=\gamma_*^\top Z+\epsilon_2,$ where ${\theta}_*\in\mathbb{R}^{d_x\times 1}$, and ${\gamma}_*\in\mathbb{R}^{d_z\times d_x}$, as also considered in Lemma (ref). Recall the true gradient in (ref) and the stochastic gradient estimator of Algorithm (ref) in (ref). Since we no longer have $X_t'$, we replace the term $X_t'$ with the predicted mean of $X_t$ given $Z_t$. Suppose that ${\gamma}_*$ is known. We specifically replace $\nabla_{{\theta_{t}}} g({\theta_{t}};X_t')=X_t'$ by $\mathbb{E}_{\mid Z_t}\left[ X_t \right] ={\gamma}_*^\top Z_t$. In such a case, indeed we have an unbiased gradient estimator:

align*[align* omitted — 528 chars of source]

where $\mathbb{E}_{t}\left[ \cdot \right] $ is the conditional expectation w.r.t the filtration defined on $\{{\gamma}_1,{\theta}_1,{\gamma}_2,{\theta}_2,\cdots,{\gamma}_t,{\theta_{t}}\}$.

In reality, ${\gamma}_*$ is unknown beforehand. Hence, we estimate ${\gamma}_*$ using some online procedure and replace $\nabla_{{\theta_{t}}} g({\theta_{t}};X_t')$ by ${\gamma_{t}}^\top Z_t$ instead of ${\gamma}_*^\top Z_t$. It leads to the following updates:

align[align omitted — 227 chars of source]

A closer inspection reveals that the updates in (ref) can diverge until ${\gamma_{t}}$ is close enough to ${\gamma}_*$. It is easy to see this fact from the following expansion of ${\theta_{t+1}}-{\theta}_*$. We have

align*[align* omitted — 504 chars of source]

where

align*[align* omitted — 177 chars of source]

However, the matrix ${\gamma_{t}}^\top\Sigma_Z{\gamma}_*$ may not be positive semi-definite, even if $\Sigma_Z$ is positive definite. Thus the negative eigenvalues associated with ${\gamma_{t}}^\top\Sigma_Z{\gamma}_*$ might cause the $\theta_t$ iterates to first diverge, before eventually converging as ${\gamma_{t}}$ gets closer to ${\gamma}_*$. We illustrate this intuition in a simple experiment in Figure (ref). To resolve this issue, we propose Algorithm (ref), where we replace $g({\theta_{t}},X_t) =X_t^\top {\theta_{t}}$ with $ Z_t^T{\gamma_{t}}{\theta_{t}}$ in (ref). With such a modification, in the corresponding decomposition for $\theta_{t+1} - \theta_*$ (see (ref)), we have ${\hat{Q}_{t}}=\left(I-{\alpha_{t+1}}{\gamma_{t}}^\top\Sigma_Z{\gamma_{t}}\right)$, where the matrix product ${\gamma_{t}}^\top\Sigma_Z{\gamma_{t}}$ is always positive semi-definite. Hence, with a properly chosen stepsize $\alpha_{t}$ we could quantify the convergence of $\theta_t$ to $\theta_*$ non-asymptotically. Nevertheless, assuming a warm-start condition on $\theta_0$, we also show the convergence of (ref), in Appendix (ref) for completeness.\\

figure[figure omitted — 416 chars of source]

Algorithm and Analysis. Based on the intuition, we present Algorithm (ref). One could interpret the algorithm as the SGD analogy of the offline 2SLS algorithm angrist1995two. It is also related to the framework of non-linear two-stage stochastic approximation algorithms doan2020finite,dalal2018finite,mokkadem2006convergence; albeit the updates of $\theta_t$ and $\gamma_t$ are coupled since both updates use $Z_t$. Furthermore, the dependency between the randomness between the two stages in the IVaR problem, makes the analysis significantly different and more challenging from the classical analysis of two-stage algorithms (see below Theorem (ref) for additional details). Finally, while Algorithm (ref) is designed for linear models, the intuition behind the method is also applicable to non-linear models (i.e., between $Z$ and $X$, and $X$ and $Y$). We focus on linear models in this work in order to derive our theoretical results. A detailed treatment of the nonlinear case (for which the analysis is significantly nontrivial) is left for future work. We make the following additional assumptions for the convergence analysis of Algorithm (ref).

algorithm[algorithm omitted — 743 chars of source]
assumptionFor some constants $C_z,C_{zy}>0$, we have the following bounds on the fourth moments: \begin{align*} \mathbb{E}_\left[ \norm{\Sigma_Z-ZZ^\top}^4 \right] \leq C_z d_z^{{\vartheta}_5}, \quad \mathbb{E}_\left[ \norm{\Sigma_{ZY}-ZY}^4 \right] \leq C_{zy}d_z^{{\vartheta}_6}, \quad \vartheta:=\max\{\vartheta_5,\vartheta_6\}.\addtocounter{equation}{1}\tag{\theequation} \end{align*}
assumptionThere exist constants $0<\mu_Z \leq \lambda_Z <\infty$ such that $\mu_Z I_{d_z}\preceq \Sigma_Z\preceq \lambda_{Z}I_{d_z} $.

The above conditions are rather mild moment conditions, similar to Assumption (ref), and could be easily verified for the linear model setting we consider.

assumption$\{{\gamma_{t}}\}_t$ is within a compact set of diameter $C_{\gamma} d_z^\varkappa$ for some constants $C_{\gamma}>0$, $\varkappa\geq 0$.

We emphasize that Assumption (ref) is only for the uncoupled sequence $\gamma_t$, which is an SGD sequence for solving a strongly-convex problem. It holds easily in various cases, for example by projecting the iterates onto any compact sets or a sufficiently large ball containing $\gamma^*$. It is also well-known that, without any projection operations, $\{{\gamma_{t}}\}_t$ sequence is almost surely bounded polyak1992acceleration under our assumptions. Finally, similar assumptions routinely appear in the analysis of SGD algorithms in various related settings; see, for example, tseng1998incremental,gurbuzbalaban2019convergence,haochen2019random,nagaraj2019sgd,ahn2020sgd,rajput2020closing.

We now present our result on the convergence of $\{{\theta_{t}}\}_t$ below in Theorem (ref) (see Appendix (ref) for the proof). In comparison to Theorem (ref) (regarding Algorithm (ref)), we highlight that Theorem (ref) provides an any-time guarantee, as the total number of iterations is not required in advance by Algorithm (ref).

theoremSuppose Assumptions (ref), (ref) (without $X'_t)$, (ref), (ref), and (ref) hold. In Algorithm (ref), for any $\iota>0$, set ${\alpha_{t}}=C_\alpha t^{-1+\iota/2}$ and ${\beta_{t}}= C_\beta t^{-1+\iota/2},$ where $C_\alpha= \min\{0.5d_z^{-4\varkappa-\vartheta/2}\lambda_Z^{-1}C_{\gamma}^{-2},0.5(\norm{{\gamma}_*}\lambda_Z)^{-2}\}$, and $C_\beta=\mu^2 d_z^{-1-2\varkappa}/128$. Then, we have \begin{align*} \mathbb{E}_\left[ \norm{{\theta_{t}}-{\theta}^*}^2 \right] =O\left(\frac{1}{t^{1-\iota}}\right). \end{align*}
remarkIn Theorem (ref), we present the step-size choices for the fastest rate of convergence. In the proof of Theorem (ref) (see Appendix (ref)), we show that convergence can be guaranteed for a range of step-sizes given by ${\alpha_{t}}=C_\alpha t^{-a}$, ${\beta_{t}}= C_\beta t^{-b},$ where $1/2<a,b<1$, $b>2-2a$ with corresponding rate being $ \mathbb{E}_{}\left[ \norm{{\theta_{t}}-{\theta}^*}^2 \right] =O(\max\{t^{-b(2-(1-\iota/2)^{-1})},t^{-a}\log(2/\iota-1)\}) $. In particular, one requires $a,b<1$ to ensure $({\alpha_{t}}-{\alpha_{t+1}})/{\alpha_{t}}=o({\alpha_{t}})$, and $({\beta_{t}}-{\beta_{t+1}})/{\beta_{t}}=o({\beta_{t}})$, as is standard in stochastic approximation literature (see, for example, chen2020statistical,polyak1992acceleration).

Proof Techniques. The major challenge towards the convergence analysis of $\{{\theta_{t}}\}_t$ lies in the interaction term ${\gamma_{t}} Z_tZ_t^\top{\gamma_{t}}{\theta_{t}}$ between $\gamma_t$ and $\theta_t$ in (ref). This multiplicative interaction term leads to an involved dependence between the noise in the stochastic gradient updates for the two stages. Such a dependence has not been considered in existing analysis of non-linear two time-scale algorithms mokkadem2006convergence,maei2009convergent,dalal2018finite,doan2020finite,xu2021sample,wang2021non,doan2022nonlinear. In addition, doan2022nonlinear considers the case when the noise sequence is not only independent of each other but also independent of iterate locations. Furthermore, they assumes (see their Assumption 3) that the condition in Assumption (ref) holds for all $\gamma$ whereas Assumption (ref) only needs to hold for ${\gamma}_*$, that is much milder. Similarly, many works (for example, Assumption 1 in wang2021non, Assumption 2 in xu2021sample and Theorem 2 in maei2009convergent) assume that the iterates of both stages are bounded in a compact set and consequently, and hence the variance of the stochastic gradients are also uniformly bounded.

In our setting, firstly, the stochastic gradient in (ref), evaluated at $({\theta_{t}},{\gamma_{t}})$ is biased:

align*[align* omitted — 568 chars of source]

Furthermore, even under Assumption (ref), the variance of the stochastic gradient is not (ref) uniformly bounded. Overcoming these issues, in addition to the aforementioned dependence between the noise in the stochastic gradient updates for the two stages, forms the major novelty in our analysis. We proceed by noting that if ${\gamma}_*$, $\Sigma_Z$, and $\Sigma_{ZY}$ were known beforehand, one conduct deterministic gradient updates, i.e., $\tilde{{\theta}}_{t+1}=\tilde{{\theta}}_t-{\alpha_{t+1}}{\gamma}_*^\top \left(\Sigma_Z{\gamma}_*\tilde{{\theta}}_t-\Sigma_{ZY}\right)$, to obtain ${\theta}_*$. By standard results on gradient descent for strongly convex functions (see, for example, nesterov2013introductory), $\{\tilde{\theta}_t\}_t$ converges exponentially fast as stated in Lemma (ref). Hence, it remains to show that the trajectory of ${\theta_{t}}$ converges to the trajectory of $\tilde{{\theta}}_t$. That is, defining the sequence ${\delta_{t}}\coloneqq {\theta_{t}}-\tilde{{\theta}}_t$, our goal is to establish the convergence rate of $\mathbb{E}_{}\left[ \lVert{\delta_{t}}\rVert_2^2 \right] $. We first provide an intermediate bound (see Lemma (ref)) and then progressively sharpen to a tighter bound (see Lemma (ref)). In doing so, it is also required to show that $\mathbb{E}_{}\left[ \norm{{\theta_{t}}}^4 \right] $ is bounded, which we prove in Lemma (ref). The proof of Lemma (ref) is non-trivial and requires carefully chosen stepsizes satisfying $\sum_{t=1}^\infty({\alpha_{t}^2}+{\alpha_{t}}\sqrt{{\beta_{t}}})<\infty$.

Numerical Experiments

Experiments for Algorithm (ref) (TOSG-IVaR). We first consider the following problem, in which $(Z, X, Y)$ is generated via

align*[align* omitted — 152 chars of source]

where $c>0$ is a scalar to control the variance of the noise vector, and $h_1$ is the first coordinate of $h$. The noise vectors (or scalar) $h, \epsilon_x, \epsilon_y$ are independent of $Z$, and we have $h\sim \mathcal{N}(\mathbf{1}_{d_x}, I_{d_x})$, $\epsilon_x \sim \mathcal{N}(0, I_{d_x})$,$\epsilon_y \sim \cN(0, 1)$. In each iteration, one tuple $(X, X', Y)$ is generated and used to update $\theta_t$ according to Algorithm (ref). We set $(d_x, d_z)\in \{(4, 8), (8, 16)\}$, $c\in \{0.1, 1.0\}$, and $\phi(s)\in \{s, s^2\}$. We repeat each setting 50 times and report the curves of $\EE[\norm{\theta_t - \theta_*}^2]$ in Figure (ref), where the expectation is computed as the average of $\norm{\theta_t - \theta_*}^2$ of all trials, and the shaded region represents the standard deviation. The first row and the second row correspond to $\phi(s) = s$ and $\phi(s) = s^2$ respectively. Here, $c=0.1$ for odd columns and $c=1.0$ for even columns. We have $(d_x, d_z) = (4, 8)$ for the first two columns and $(d_x, d_z) = (8, 16)$ for the last two columns. Empirically, we can observe that our Algorithm (ref) performs well across all different settings.

figure[figure omitted — 1,199 chars of source]

\noindentExperiments for Algorithm (ref) (OTSG-IVaR). Next, we compare our Algorithm (ref) as well as its variant and Algorithm 1 in della2023online. We write “OTSG-IVaR”, “CSO -- Eq. (11)” and “[DVB23]” to represent Algorithm (ref), Algorithm (ref) with the updates replaced by (ref) and Algorithm 1 in della2023online (see Appendix (ref)). We follow simulation settings similar to della2023online:

align*[align* omitted — 234 chars of source]

where $\epsilon_1$ is the first coordinate of $\epsilon$, $\theta_*\in\mathbb{R}^{d_x}$ is a unit vector chosen uniformly randomly, and $\gamma_*\in\mathbb{R}^{d_z\times d_x}$ where $\gamma_{ij}=0$ for $i\neq j$, and $\gamma_{ij}=1$ for $i= j$, $i=1,2,\cdots,d_x$, and $j=1,2,\cdots,d_z$. Here $\rho$ controls the level of endogeneity in the model. We compare the performance of Algorithm (ref) with (ref), and O2SLS della2023online for $\rho=1,4$, and $\sigma_\epsilon=0.5,1$. By varying $\sigma_\epsilon$ we control the correlation between $X$ and $Z$. We consider two settings $(d_x,d_z)=(1,1)$, and $(d_x,d_z)=(8,16)$. As performance metric, in Figure (ref) we plot $\mathbb{E}_{}\left[ \norm{{\theta_{t}}-{\theta}_*}^2 \right] $ where the $\mathbb{E}_{}\left[ \cdot \right] $ is approximated by averaging over 50 trials, and both axes are in $\log$ scale (base 10). We also show, in Figure (ref), the convergence of the test Mean Squared Error (MSE) evaluated over 400 test samples to the best possible Test MSE where ${\theta}_*$ and ${\gamma}_*$ are known beforehand. For Figures (ref) and (ref), the first row and second row corresponds to $(d_x, d_z) = (1, 1)$ and $(d_x, d_z) = (8, 16)$ respectively, and $\sigma_{\epsilon} = 0.5$ in odd columns and $\sigma_{\epsilon} = 1.0$ in even columns. We have $\rho = 1.0$ for the first two columns and $\rho = 4.0$ for the last two columns. We can observe that O2SLS has much larger variance in different settings, while our algorithms perform consistently well in all settings.

figure[figure omitted — 1,459 chars of source]
figure[figure omitted — 1,255 chars of source]

Conclusion

We presented streaming algorithms for least-squares IVaR based on directly solving the associated conditional stochastic optimization formulation in (ref). Our algorithms have several benefits, including avoidance of mini-batches and matrix inverses. We show that the expected rates of convergences for the proposed algorithms are of order $\mathcal{O}(\log T/T)$ and $\mathcal{O}(1/T^{1-\iota})$, for any $\iota>0$, under the availability of two-sample and one-sample oracles, respectively. As future work, it is interesting to develop streaming inferential methods for IVaR. Leveraging related works for the vanilla SGD polyak1992acceleration, anastasiou2019normal, shao2022berry, chen2020statistical, zhu2023online to the setting of Algorithms (ref) and (ref), provides a concrete direction to establish Central Limit Theorems and develop limiting covariance estimation procedures.

appendix