Stochastic Optimization Algorithms for Instrumental Variable Regression with Streaming Data

Authors: xuxing  chen
abhishek  roy
yifan  hu
krishnakumar  balasubramanian
Venue: NeurIPS_cc_2024_Conference
Published: 2024-04-24 13:41:41

Abstract: We 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 thereby providing 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-\epsilon})$ for any $\epsilon>0$, respectively under the availability of two-sample and one-sample oracles respectively. Importantly, under the availability of the two-sample oracle, the aforementioned rate is actually agnostic to the relationship between confounder and the instrumental variable demonstrating the flexibility of the proposed approach in alleviating the need for explicit model assumptions required in recent works based on reformulating the problem as min-max optimization problems. Experimental validation is provided to demonstrate the advantages of the proposed algorithms over classical approaches like the 2SLS method.

Main paper: 
 introduction: 
 Instrumental variable analysis is widely used in fields like econometrics, health care [TM17] , social science [Bol12] , 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:Y = g θ * (X) + ϵ 1 with X = h γ * (Z) + ϵ 2 (1)where X ∈ R dx and ϵ 1 are correlated and ϵ 2 is a centered unobserved noise (independent of Z ∈ R dz ), leading to confounding in the model between X and Y ∈ R. Here ϵ 1 and ϵ 2 are dependent, and θ * and γ * are true parameters for the respective function g and h. Our goal is to design efficient algorithms that recovers θ * from the data.♯ XC and AR contributed equally to this work.Traditionally, IVaR algorithms are based on two-stage estimation procedures, where we first regress Z and X to obtain an estimator X, and then regress X and Y , with the essence that 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 [HH05, DFFR11, HLLBT17], with the parametric two-stage least-squares (2SLS) procedure being the most canonical one [AI95] .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, [AP09] call the first stage regression as "forbidden regression", due to the concerns in estimating a nuisance parameter.Considering the squared loss function, [MMLR20] formulate the IVaR problem as a conditional stochastic optimization problem [HZCH20] :EQUATIONHowever, [MMLR20] did not solve problem (2) efficiently, and resort to reformulating (2) further as a minimax optimization problem. Indeed, they mention explicitly in their work that "it remains cumbersome to solve (2) 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 [DHP + 17], [MMLR20] propose to use reproducing kernel Hilbert space (RKHS) to handle the maximization over continuous functional space. See also [LS18, BKS19, DLMS20, LCY + 20, BKM + 23] 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 [MMLR20] that problem (2) is cumbersome to solve, we design and analyze efficient streaming algorithms to directly solve the conditional stochastic optimization problem in (2). Direct application of methods from [HZCH20] for solving (2) is possible, yet their approach utilizes nested sampling, i.e., for each sample of Z, [HZCH20] generate a batch of samples of X from P(X|Z), to reduce the bias in estimating the composition of nonlinear loss function with conditional expectations. Thus their methods are not suitable for the streaming setting that we are interested in. Considering (2), we first parameterize the function class G := {g(θ; X) | θ ∈ R d θ }. Now, defining F (g) := F (θ), we observe that the gradient ∇F (θ) admits the following formEQUATIONwhich implies that one does not need the nested sampling technique to reduce the bias. However, the presence of product of two conditional expectations E X|Z [g(θ; 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 P(Z|X).Contributions. We make the following contributions in this work.• 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 (3). Based on this crucial observation, we propose the Two-Sample One-stage Stochastic Gradient IVaR (TOSG-IVaR) method (Algorithm 1) 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 O(log T /T ) for the proposed method, where T is the overall number of iterations; see Theorem 1. • One-sample oracles: In the case when we do not have the aforementioned two-sample oracle, we estimate the stochastic gradient in (3) 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 2). Assuming further that the X depends linearly on the instrument Z, we establish a rate of convergence of order O(1/T 1-ι ), for any ι > 0; see Theorem 2.1.1 Literature Review IVaR analysis. Instrumental variable analysis has a long history, starting from the early works by [Wri28] and [Rei45] . Several works considered the aforementioned two-stage procedure for IVaR; a summary could be found in the work by [AP09] . Nonparametric approaches based on wavelets, splines, reproducing kernels and deep neural networks could be found, for example, in the works by [HLLBT17, SSG19, BKS19, MMLR20 Sample average approximation algorithms in this context are considered in the works of [EN13] and [HCH20] . Optimal iterative stochastic optimization algorithms for the case of T = 2 were by derived by [GRW20] . For the general T ≥ 1 case, [WFL17] provided sub-optimal rates, whereas [BGN22] derived optimal rates; see also [ZX21] and [CSY21] for related works under stronger assumptions, and [Rus21] for similar asymptotic results. While the above works required certain independence assumptions regarding the randomness across the different compositions, [HZCH20, HWX + 24] 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 (2) falls in. Compared to prior works, for e.g., [GRW20] and [BGN22] , in order to handle the dependency between the levels, [HZCH20] 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 (2) being a conditional stochastic optimization problem, mini-batches are not required due the additional favorable quadratic structure available in IVaR.Streaming IVaR. [VSH + 16, DVB24] analyzed streaming versions of 2SLS in the online 1 and adversarial settings. Focusing on linear models, [VSH + 16] provide preliminary asymptotic analysis assuming access to efficient no-regret learners, while [DVB24] 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 [DVB24] ; see Sections A and B for details. [CLL + 23] 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 (2). By [MMLR20, Theorem 4], the optimal solution of (2) gives the true underlying causal relationship under the following assumption.Assumption 2.1. (Identifiability Assumption)• The conditional distribution P Z|X is continuous in Z for any value of X. Update θ t θ t+1 = θ t -α t+1 (g(θ t ; X t ) -Y t )∇ θ g(θ t ; X ′ t ). 4: end for Output: θ T . • The function class G := {g(θ; X) | θ ∈ R d θ }is correctly specified, i.e., it includes the true underlying relationship between X and Y .Notice that both assumptions are standard in the IVaR literature [NP03, CP12, MMLR20], and makes the objective in (2) is the meaningful for IVaR. However, [MMLR20] resort to reformulating the objective function in (2) as a minimax optimization problem as described in Section 1. 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 (2), leveraging the structure provided by the quadratic loss. Given the gradient formulation in (3), a natural way to build unbiased gradient estimator is to generate X and X ′ , two independent samples of X from the conditional distributions P X|Z , for a given realization of Z and generate one sample of Y from the conditional distribution P Y |X . Then, an unbiased gradient estimator isEQUATIONThis 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 1. In particular, the algorithm never requires estimating (or modeling) the relationship between X and Z as needed in the two-stage procedure [AP09] and the minimax formulation based procedures [MMLR20, LS18, BKS19, DLMS20, LCY + 20, BKM + 23]. 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.Assumption 2.2. The tuple (Z t , X t , X ′ t , Y t ) is independent and identically distributed, across t.To demonstrate the convergence rate of Algorithm 1, we first consider the case when g is a linear function, i.e., g(θ; X) = X ⊤ θ. We make the following assumptions.Assumption 2.3. Suppose there exists µ > 0 such thatE Z E X|Z [X] • E X|Z [X] ⊤ ⪰ µI. Assumption 2.4. Let (ϑ 1 , ϑ 2 , ϑ 3 , ϑ 4 ) ∈ R 4 + .For any Z, X ′ and X i.i.d. generated from P Z|X , andEQUATIONEQUATIONEQUATIONEQUATIONwhere ∥ • ∥ 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 2.3 and 2.4 are naturally satisfied even under non-linear modeling assumption on (1). We defer its proof to Section D. := C y d ϑ2 x + C yx d ϑ4 z , set α t ≡ α = log T µT ≤ µ µ 2 +3σ 2 1. Then, we haveE ∥θ T -θ * ∥ 2 ≤ E ∥θ 0 -θ * ∥ 2 T + 3∥θ * ∥ 2 (σ 2 1 + σ 2 2 ) log T µ 2 T .Proof techniques. In the analysis of Theorem 1, the following decomposition (see (18) for the derivation) plays a crucial role:θ t+1 -θ * = A t + α t+1 B t , A t = θ t -α t+1 E Z E X|Z [X] • E X|Z [X] ⊤ θ t + α t+1 E Z E Y |Z [Y ] • E X|Z [X] -θ * , B t = -X ′ t X ⊤ t -E Z E X|Z [X] • E X|Z [X] ⊤ θ t + Y t X ′ t -E Z E Y |Z [Y ] • E X|Z [X] ,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 [Lan20] and the expected smoothness condition [KR20] . In the IVaR setup, such standard assumptions do not hold as θ 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 [MMLR20] , notice that we use two samples of X from the conditional distribution P X|Z and achieve an O(1/T ) last iterate convergence rate to the global optimal solution, which is the true underlying causal relationship under Assumption 2.1. In comparison, [MMLR20] 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 [CF21] ), potentially using the Multilevel Monte Carlo (MLMC) based stochastic gradient estimators. While [HCH21] , 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(θ; 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. Assumption 2.5. Let the following assumptions hold:• Function F (θ) is ℓ-smooth.• The iterates {θ t } T +1 t=1 generated by Algorithm 1 are in a compact set A. • The random objects X|Z and Y |Z have bounded variance for any Z, i.e., there exist σ > 0 such thatE ∥X -E [X | Z] ∥ 2 | Z ≤ σ 2 , E ∥Y -E [Y | Z] ∥ 2 | Z ≤ σ 2 .Proposition 1. Suppose Assumptions 2.1, 2.2, and 2.5 hold. Choosingα t ≡ α = O 1 √ T , for Algorithm 1 we have min 1≤t≤T E ∥∇F (θ t )∥ 2 = O 1 √ T .The proof of the proposition is immediate. Note that under Assumption 2.5, we can deduce that the unbiased gradient estimator v(θ) = (g(θ; X) -Y )∇ θ g(θ; X ′ ) has a bounded variance sinceVar(v(θ)) =Var(g(θ; X) -Y )Var(∇ θ g(θ; X ′ )) + Var(g(θ; X) -Y ) (E [∇ θ g(θ; X ′ )]) 2 + Var(∇ θ g(θ; X ′ )) (E [g(θ; X) -Y ]) 2 ≤ σ 2 v ,where the variance and expectation are taken conditioning on Z and θ, and σ v > 0 is a constant that only depends on σ, function g and the compact set A in Assumption 2.5. Then one can directly follow the analysis of non-convex stochastic optimization (see, for example, [GL13, Theorem 2.1]) to obtain Proposition 1. Relaxing the Assumption 2.5 (typically made in the stochastic optimization literature) with more natural assumptions on the statistical model and obtaining a result as in Theorem 1 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 (4) makes the stochastic gradient estimator biased.Intuition. Consider the case of linear models, i.e., Y = θ ⊤ * X + ϵ 1 with X = γ ⊤ * Z + ϵ 2 , where θ * ∈ R dx×1 , and γ * ∈ R dz×dx , as also considered in Lemma 1. Recall the true gradient in (3) and the stochastic gradient estimator of Algorithm 1 in (4). 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 γ * is known. We specifically replace∇ θt g(θ t ; X ′ t ) = X ′ t by E |Zt [X t ] = γ ⊤ * Z t .In such a case, indeed we have an unbiased gradient estimator:E t γ ⊤ * Z t (X ⊤ t θ t -Y t ) = E t E |Zt [X t ] (E |Zt [X t ] ⊤ θ t -E |Zt [Y t ]) =E t γ ⊤ * Z t Z ⊤ t γ * (θ t -θ * ) = γ ⊤ * Σ Z γ * (θ t -θ * ) = ∇ θ F (θ t ), where E t [•] is the conditional expectation w.r.t the filtration defined on {γ 1 , θ 1 , γ 2 , θ 2 , • • • , γ t , θ t }.In reality, γ * is unknown beforehand. Hence, we estimate γ * using some online procedure and replace∇ θt g(θ t ; X ′ t ) by γ t ⊤ Z t instead of γ ⊤ * Z t .It leads to the following updates:EQUATIONA closer inspection reveals that the updates in (11) can diverge until γ t is close enough to γ * . It is easy to see this fact from the following expansion of θ t+1 -θ * . We haveEQUATIONwhereξ Zt = Σ Z -Z t Z ⊤ t , ξ ZtYt = Σ ZY -Z t Y t , Q t := I -α t+1 γ t ⊤ Σ Z γ * .However, the matrix γ t ⊤ Σ Z γ * may not be positive semi-definite, even if Σ Z is positive definite. Thus the negative eigenvalues associated with γ t ⊤ Σ Z γ * might cause the θ t iterates to first diverge, before .EQUATIONEQUATION4: end for eventually converging as γ t gets closer to γ * . We illustrate this intuition in a simple experiment in Figure 1 . To resolve this issue, we propose Algorithm 2, where we replace g(θ t , X t ) = X ⊤ t θ t with Z T t γ t θ t in (11). With such a modification, in the corresponding decomposition for θ t+1 -θ * (see (40)), we have Q t = I -α t+1 γ t ⊤ Σ Z γ t , where the matrix product γ t ⊤ Σ Z γ t is always positive semi-definite. Hence, with a properly chosen stepsize α t we could quantify the convergence of θ t to θ * non-asymptotically. Nevertheless, assuming a warm-start condition on θ 0 , we also show the convergence of (11), in Appendix E.3 for completeness.Algorithm and Analysis. Based on the intuition, we present Algorithm 2. One could interpret the algorithm as the SGD analogy of the offline 2SLS algorithm [AI95] . It is also related to the framework of non-linear two-stage stochastic approximation algorithms [DR20, DTSM18, MP06]; albeit the updates of θ t and γ 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 2 for additional details). Finally, while Algorithm 2 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 2. Assumption 3.1. For some constants C z , C zy > 0, we have the following bounds on the fourth moments:EQUATION) Assumption 3.2. There exist constants 0 < µ Z ≤ λ Z < ∞ such that µ Z I dz ⪯ Σ Z ⪯ λ Z I dz .The above conditions are rather mild moment conditions, similar to Assumption 2.4, and could be easily verified for the linear model setting we consider. Assumption 3.3. {γ t } t is within a compact set of diameter C γ d κ z for some constants C γ > 0, κ ≥ 0.We emphasize that Assumption 3.3 is only for the uncoupled sequence γ 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 γ * . It is also well-known that, without any projection operations, {γ t } t sequence is almost surely bounded [PJ92] under our assumptions. Finally, similar assumptions routinely appear in the analysis of SGD algorithms in various related settings; see, for example, [Tse98, GOP19, HS19, NJN19, AYS20, RGP20].We now present our result on the convergence of {θ t } t below in Theorem 2 (see Appendix E.1 for the proof). In comparison to Theorem 1 (regarding Algorithm 1), we highlight that Theorem 2 provides an any-time guarantee, as the total number of iterations is not required in advance by Algorithm 2. Theorem 2. Suppose Assumptions 2.3, 2.2 (without X ′ t ), 3.1, 3.3, and 3.2 hold. In Algorithm 2, for any ι > 0, set α t = C α t -1+ι/2 and β t = C β t -1+ι/2 , where C α = min{0.5d-4κ-ϑ/2 z λ -1 Z C -2 γ , 0.5(∥γ * ∥λ Z ) -2 }, and C β = µ 2 d -1-2κ z /128. Then, we haveE ∥θ t -θ * ∥ 2 = O 1 t 1-ι .Remark 1. In Theorem 2, we present the step-size choices for the fastest rate of convergence. In the proof of Theorem 2 (see Appendix E.1), we show that convergence can be guaranteed for a range of step-sizes given by αt = C α t -a , β t = C β t -b , where 1/2 < a, b < 1, b > 2 -2a with corresponding rate being E ∥θ t -θ * ∥ 2 = O(max{t -b(2-(1-ι/2) -1 ) , t -a log(2/ι -1)}). In particular, one requires a, b < 1 to ensure (α t -α t+1 )/α t = o(α t ), and (β t -β t+1 )/β t = o(β t ),as is standard in stochastic approximation literature (see, for example, [CLTZ20, PJ92] ).Proof Techniques. The major challenge towards the convergence analysis of {θ t } t lies in the interaction term γ t Z t Z ⊤ t γ t θ t between γ t and θ t in (13). 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 [MP06, MSB + 09, DTSM18, DR20, XL21, WZZ21, Doa22]. In addition, [Doa22] 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 2.3 holds for all γ whereas Assumption 2.3 only needs to hold for γ * , that is much milder. Similarly, many works (for example, Assumption 1 in [WZZ21] , Assumption 2 in [XL21] and Theorem 2 in [MSB + 09]) 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 (13), evaluated at (θ t , γ t ) is biased:E t,Zt γ t ⊤ Z t (Z ⊤ t γ t θ t -Y t ) = E t,Zt γ t ⊤ Z t (Z ⊤ t γ t θ t -Z ⊤ t γ * θ * ) = E t γ t ⊤ Σ Z (γ t θ t -γ * θ * ) =γ t ⊤ Σ Z γ t (θ t -θ * ) + γ t ⊤ Σ Z (γ t -γ * )θ * ̸ = γ ⊤ * Σ Z γ * (θ t -θ * ) = ∇ θ F (θ t ). Furthermore, even under Assumption 3.3, the variance of the stochastic gradient is not (13) 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 γ * , Σ Z , and Σ ZY were known beforehand, one conduct deterministic gradient updates, i.e., θ t+1 = θ t -α t+1 γ ⊤ * Σ Z γ * θ t -Σ ZY , to obtain θ * . By standard results on gradient descent for strongly convex functions (see, for example, [Nes13] ), { θ t } t converges exponentially fast as stated in Lemma 4. Hence, it remains to show that the trajectory of θ t converges to the trajectory of θ t . That is, defining the sequence δ t := θ t -θ t , our goal is to establish the convergence rate of E ∥δ t ∥ 2 2 . We first provide an intermediate bound (see Lemma 6) and then progressively sharpen to a tighter bound (see Lemma 7). In doing so, it is also required to show that E ∥θ t ∥ 4 is bounded, which we prove in Lemma 5. The proof of Lemma 5 is non-trivial and requires carefully chosen stepsizes satisfying∞ t=1 (α 2 t + α t √ β t ) < ∞. 
numerical experiments: 
 Experiments for Algorithm 1 (TOSG-IVaR). We first consider the following problem, in which (Z, X, Y ) is generated via 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, ϵ x , ϵ y are independent of Z, and we have h ∼ N (1 dx , I dx ), ϵ x ∼ N (0, I dx ),ϵ y ∼ N (0, 1). In each iteration, one tuple (X, X ′ , Y ) is generated and used to update θ t according to Algorithm 1. We set (d x , d z ) ∈ {(4, 8), (8, 16)}, c ∈ {0.1, 1.0}, and ϕ(s) ∈ {s, s 2 }. We repeat each setting 50 times and report the curves of E[∥θ t -θ * ∥ 2 ] in Figure 2 , where the expectation is computed as the average of ∥θ t -θ * ∥ 2 of all trials, and the shaded region represents the standard deviation. The first row and the second row correspond to ϕ(s) = s and ϕ(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 1 performs well across all different settings.Z ∼ N (0, I dz ), X = ϕ(γ ⊤ * Z) + c • (h + ϵ x ), Y = θ ⊤ * X + c • (h 1 + ϵ y ),Experiments for Algorithm 2 (OTSG-IVaR). Next, we compare our Algorithm 2 as well as its variant and Algorithm 1 in [DVB24] . We write "OTSG-IVaR", "CSO -Eq. ( 11)" and " [DVB23] " to represent Algorithm 2, Algorithm 2 with the updates replaced by (11) and Algorithm 1 in [DVB24] (see Appendix A). We follow simulation settings similar to [DVB24] :EQUATIONwhere ϵ 1 is the first coordinate of ϵ, θ * ∈ R dx is a unit vector chosen uniformly randomly, and γ * ∈ R dz×dx where γ ij = 0 for i ̸ = j, andγ ij = 1 for i = j, i = 1, 2, • • • , d x , and j = 1, 2, • • • , d z .Here ρ controls the level of endogeneity in the model. We compare the performance of Algorithm 2 with (11), and O2SLS [DVB24] for ρ = 1, 4, and σ ϵ = 0.5, 1. By varying σ ϵ 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 3 we plot E ∥θ t -θ * ∥ 2 where the E [•] is approximated by averaging over 50 trials, and both axes are in log scale (base 10). We also show, in Figure 4 , the convergence of the test Mean Squared Error (MSE) evaluated over 400 test samples to the best possible test MSE where θ * and γ * are known beforehand. For Figures 3 and 4 , the first row and second row corresponds to (d x , d z ) = (1, 1) and (d x , d z ) = (8, 16) respectively, and σ ϵ = 0.5 in odd columns and σ ϵ = 1.0 in even columns. We have ρ = 1.0 for the first two columns and ρ = 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. We further conduct experiments on real-world datasets provided in [AE96] and [Rya12] . Due to space limit, we include the numerical results in Section C.3 of the Appendix.t ∥ 2 < 1, ∥Σ ZY ∥ 2 = O(1), Cauchy-Schwarz inequality, and Lemma 3, we get,EQUATIONCombining (31), (32), (33), (34), (35), and Lemma 5, we haveEQUATIONIn the above, the third inequality follows by choosing β t ≤ µ 2 /(64d z ), and α t+1 √ d z β t < 1. Then, from (37), we haveE ∥δ t ∥ 2 2 = O β t d 1+2κ z + α t+1 d 4κ+ϑ/2 z + d z β t .Coming back to the proof of Theorem 2, observe that, we can sharpen the bound in (35) using Lemma 6 which allows us to avoid the use of Young's inequality. This leads to the following improved version of the recursion in (37) using which we can improve the term √ d z β t in (29) as follows:E ∥δ t+1 ∥ 2 2 ≲(1 -7µα t+1 /8)E ∥δ t ∥ 2 + α t+1 O β t d 1+2κ z + α t+1 d 4κ+ϑ/2 z + α t+1 β t d 1/2+2κ+ϑ/4 z + (β t d z ) 3/4 =O β t d 1+2κ z + α t+1 d 4κ+ϑ/2 z + α t+1 β t d 1/2+2κ+ϑ/4 z + (β t d z ) 3/4 .In fact, this trick can be used repeatedly to sharpen the bound even further as shown in Lemma 7.Lemma 7 (Final improved bound on E[∥δ t ∥ 2 2 ]). Let the conditions in Theorem 2 be true. Then using Lemma 6, we have,E ∥δ t+1 ∥ 2 2 ≲O (d z β t ) 1-2 -r-1 + r i=0 α 2 -i t+1 β t 1-2 -i d 1+(4κ+ϑ/2-1)2 -i z + β t (1 + α 2 -i t+1 )d 1+2 1-i κ z ,where r is any non-negative integer.Proof of Lemma 7. If we haveE ∥δ t ∥ 2 = O α t+1 d 4κ+ϑ/2 z + β t d 1+2κ z + d z β t ,then from (35), we have,EQUATIONThen, similar to (36), we have,E ∥δ t+1 ∥ 2 2 ≲(1 -α t+1 µ)E ∥δ t ∥ 2 + 4α 2 t+1 β t d 1+2κ z + 4Cα 2 t+1 d 4κ+ϑ/2 z + 2α t+1 µE ∥δ t ∥ 2 /16 + 4d 1+2κ z β t /µ + α t+1 β t d 1/2+2κ+ϑ/4 z + β t d 1+κ z + (d z β t ) 3/4 ≲(1 -7µα t+1 /8)E ∥δ t ∥ 2 + α t+1 O (d z β t ) 3/4 + 1 i=0 α 2 -i t+1 β t 1-2 -i d 1+(4κ+ϑ/2-1)2 -i z + β t (1 + α 2 -i t+1 )d 1+2 1-i κ z =O (d z β t ) 3/4 + 1 i=0 α 2 -i t+1 β t 1-2 -i d 1+(4κ+ϑ/2-1)2 -i z + β t (1 + α 2 -i t+1 )d 1+2 1-i κ z .Now if we repeat this step r number of times (where r is to be set later), by progressive sharpening we get the following bound.E ∥δ t+1 ∥ 2 2 ≲O (d z β t ) 1-2 -r-1 + r i=0 α 2 -i t+1 β t 1-2 -i d 1+(4κ+ϑ/2-1)2 -i z + β t (1 + α 2 -i t+1 )d 1+2 1-i κ z .Coming back to the proof of Theorem 2, we have that by combining Lemma 4, and Lemma 7,EQUATIONNow, in (39), for some arbitrarily small number ι > 0, choosingα t = min(0.5d -4κ-ϑ/2 z λ -1 Z C -2 γ , 0.5(∥γ * ∥ 2 λ Z ) -2 )t -1+ι/2 , β t = µ 2 d -1-2κ z t -1+ι/2 /128,and setting r = ⌈log 2 ((ι/2) -1 -1) -1⌉ we get,E ∥θ t -θ * ∥ 2 = O max t -1+ι , t -1+ι/2 log((ι/2) -1 -1) . 
conclusion: 
 We presented streaming algorithms for least-squares IVaR based on directly solving the associated conditional stochastic optimization formulation in (2). 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 O(log T /T ) and O(1/T 1-ι ), for any ι > 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 [PJ92, ABE19, SZ22, CLTZ20, ZCW23] to the setting of Algorithms 1 and 2, provides a concrete direction to establish Central Limit Theorems and develop limiting covariance estimation procedures. + d z ) + d z d x + d x O(d 3 x + td 2 x + td x d z ) O2SLS [DVB24, v3] (Sec. A) d 2 x + d 2 z + d z d x + d x O(d 2 x + d 2 z + d z d x ) TOSG-IVaR (our Alg 1) d x O(d x )OTSG-IVaR (our Alg 2)d x d z + d x O(min(d 2 x , d 2 z ) + d z d x )C Experimental Details 
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: Sections 2 and 3. Our work is mainly about theory and algorithm design, and 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? 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. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: All experiments in Section 4 were conducted on a computer with an 11th Intel(R) Core(TM) i7-11370H CPU. The time and space required to run our experiments are negligible and we anticipate they can be conducted in almost all computers.C.2 Experimental Details for Figure 1 In Figure 1 , we show an example where the updates (11) may diverge first before converging eventually and finite time performance can be much worse compared to Algorithm 2. For this experiment, we choose the model presented in (16) withd x = d z = 1, θ * = 1, γ * = -1, ρ = 4, and σ ϵ = 1. When initialized at γ 0 = 10, and θ 0 = 0, the updates in (11) keeps diverging rapidly at first whereas Algorithm 2 is much more stable. So, by the end of 100, 000 iterations, while Algorithm 2 achieves an error of ≈ 10 -5 , (11) achieves ≈ 10 4 that is worse than it was at initialization because (11) has not recovered from the initial divergence phase yet. However, once (11) starts converging, the convergence rate of (11) is similar to Algorithm 2 as one can see from Figure 1 (also see our discussion on the convergence of (11) in Section E.3).In this section we provide experimental results on real-world datasets provided in [AE96] and [Rya12] .The results are included in Figures 5 and 6 respectively. Following the convention in Section 4, we write "OTSG-IVaR", "CSO -Eq. ( 11)" and " [DVB23] " to represent Algorithm 2, Algorithm 2 with the updates replaced by (11) and Algorithm 1 in [DVB24] .The outcome Y is number of working weeks divided by 52, the regressor X is 1(number of children is greater than 2), and the instrumental variable Z is 1(first two siblings are of same sex), where 1(•) is the indicator function. At each time we randomly sample a data-point from this dataset without replacement. Since the "true" parameter θ * is unknown in real data, we use the offline model parameter estimate as our ground truth, following [CLL + 23]. We include the results in Figure 5 , in which we observe that CSO performs similar to [DVB23] whereas OTSG-IVaR performs much better in terms of convergence speed. Moreover, the estimation errors of CSO and [DVB23] plateau after ≈ 10,000 iterations whereas OTSG-IVaR keeps on improving the estimate over the observed horizon.Response variable Y is log(shipped), and the predictor X is log(price). There are 4 instrumental variables given by wage in dollars per hour for skilled manufacturing workers, electricity price, coal price, and gas price. Unlike the previous dataset, we have only 483 data samples here. So, to mimic an i.i.d. data stream, we divide our training into multiple epochs of length equal to the number of training data samples, and over each epoch, at each iteration, we sample one data point uniformly randomly without replacement from the training data to generate the data-stream. Both OTSG-IVaR and CSO perform much better than [DVB23] in terms of convergence speed. Figure is a magnified view of iterations > 50 to highlight the performance difference between various algorithms. Proof. We aim to find the optimal θ * . According to (2), we knowEQUATIONThe updates in Algorithm 1 can be written asEQUATIONNow we analyze the convergence and variance separately. For the convergence part, we haveEQUATIONFor the variance part we haveE ∥X ′ t X ⊤ t -E Z E X|Z [X] • E X|Z [X] ⊤ ∥ 2 =E ∥X ′ t X ⊤ t -E X|Zt [X] • E X|Zt [X] ⊤ + E X|Zt [X] • E X|Zt [X] ⊤ -E Z E X|Z [X] • E X|Z [X] ⊤ ∥ 2 ≤2E ∥X ′ t X ⊤ t -E X|Zt [X] • E X|Zt [X] ⊤ ∥ 2 + ∥E X|Zt [X] • E X|Zt [X] ⊤ -E Z E X|Z [X] • E X|Z [X] ⊤ ∥ 2 ≤2C x d ϑ1 x + 2C xx d ϑ3 z =: σ 2 1 (20)Similarly, we haveEQUATION)Now we know from ( 18), ( 19), (20), and (21) thatEQUATION)whereA t = I -α t+1 E Z E X|Z [X] • E X|Z [X] ⊤ (θ t -θ * ) B t = -X ′ t X ⊤ t -E Z E X|Z [X] • E X|Z [X] ⊤ θ t + Y t X ′ t -E Z E Y |Z [Y ] • E X|Z [X] .This impliesEQUATION)where the first inequality uses Cauchy-Schwarz inequality, the definition of σ 1 , σ 2 and Assumption 2.4. Choosing α t+1 such that((1 -α t+1 µ) 2 + 3α 2 t+1 σ 2 1 ) ≤ 1 -α t+1 µ ⇔ α ≤ µ µ 2 + 3σ 2 1and taking expectation on both sides of (23), we haveE ∥θ t+1 -θ * ∥ 2 ≤ (1 -α t+1 µ)E ∥θ t -θ * ∥ 2 + 3α 2 t+1 σ 2 1 ∥θ * ∥ 2 + 3α 2 t+1 σ 2 2 ∥θ * ∥ 2 .Now, we use the following result.Lemma 2. Suppose we have three sequencest } ∞ t=0 , {b t } ∞ t=0 , {r t } ∞ t=0 satisfying a t+1 ≤ r t a t + b t , r t > 0 (24)for any t ≥ 0. Define R t+1 = t i=0 r i , we havea t+1 ≤ R t+1 a 0 + t i=0 R t+1 b i R i+1 .By Lemma 2, we knowE ∥θ t+1 -θ * ∥ 2 ≤ t i=0 (1 -α i µ)E ∥θ 0 -θ * ∥ 2 + (3σ 2 1 ∥θ * ∥ 2 + 3σ 2 2 ∥θ * ∥ 2 ) t i=0 α 2 i t j=i+1(1 -α j µ).i = α, we know E ∥θ t -θ * ∥ 2 ≤(1 -αµ) t E ∥θ 0 -θ * ∥ 2 + α 2 t i=0 (1 -αµ) i (3σ 2 1 ∥θ * ∥ 2 + 3σ 2 2 ∥θ * ∥ 2 ) ≤e -tαµ E ∥θ 0 -θ * ∥ 2 + α µ (3σ 2 1 ∥θ * ∥ 2 + 3σ 2 2 ∥θ * ∥ 2 ) Choosing α, T such that α = log T µT ≤ µ µ 2 +3σ 2 1 , we know E ∥θ T -θ * ∥ 2 ≤ E ∥θ 0 -θ * ∥ 2 T + 3∥θ * ∥ 2 (σ 2 1 + σ 2 2 ) log T µ 2 T .Proof. We notice from (24) that for any 0 ≤ i ≤ t, we havea i+1 R i+1 ≤ a i R i + b i R i+1 .Taking summation on both sides, we havea t+1 R t+1 ≤ a 0 R 0 + t i=0 b i R i+1which completes the proof by multiplying R t+1 on both sides.Proof. We first notice that Assumption 2.3 holds since 5) and (6), we haveE Z E X|Z [X] • E X|Z [X] ⊤ = E Z ϕ(γ ⊤ * Z) • ϕ(γ ⊤ * Z) ⊤ ⪰ µI. For (EQUATIONandE ∥Y X ′ -E Y |Z [Y ]E X|Z [X]∥ 2 =E ∥X ′ X ⊤ θ * + ϵ 1 X ′ -E X|Z [X]E X|Z [X] ⊤ θ * ∥ 2 ≤2E ∥X ′ X ⊤ θ * -E X|Z [X]E X|Z [X] ⊤ θ * ∥ 2 + 2E ϵ 2 1 ∥ϕ(γ ⊤ * Z) + ϵ ′ 2 ∥ 2 =O(∥θ * ∥ 2 σ 2 ϵ2 d 2 x + σ 2 ϵ1 d x + σ 2 ϵ1,ϵ2 d x ), where the first inequality uses Cauchy-Schwarz inequality, and the second equality uses (9), ( 10) and (25). For (7) we haveE ∥E X|Z [X] • E X|Z [X] ⊤ -E Z E X|Z [X] • E X|Z [X] ⊤ ∥ 2 =E ∥ϕ(γ ⊤ * Z)ϕ(γ ⊤ * Z) ⊤ -E ϕ(γ ⊤ * Z)ϕ(γ ⊤ * Z) ⊤ ∥ 2 = O(d z )where the last equality uses (10). Using the above conclusion in (8), we haveE ∥E Y |Z [Y ] • E X|Z [X] -E Z E Y |Z [Y ] • E X|Z [X] ∥ 2 =E ∥E X|Z [X] • E X|Z [X] ⊤ θ * -E Z E X|Z [X] • E X|Z [X] ⊤ θ * ∥ 2 = O(∥θ * ∥ 2 d z ).Proof of Theorem 2 . Recall that ξ Zt , and ξ ZtYt are the i.i.d. noise sequencesξ Zt = Σ Z -Z t Z ⊤ t , ξ ZtYt = Σ ZY -Z t Y t .Note γ * , and θ * can be written as γ * = Σ -1 Z Σ ZX ∈ R dz×dx , and θ * = γ ⊤ * Σ z γ * -1 γ ⊤ * Σ ZY ∈ R dx which we are going to use throughout the proof.To quantify the bias, we use the following bound onE ∥γ t -γ * ∥ k 2 , k = 1, 2, 4, proved in Lemma 3.2 of [CLTZ20].Lemma 3. Suppose Assumption 2.2, and Assumption 3.2 hold. Then we haveEQUATIONWe proceed by noting that if γ * , Σ Z , and Σ ZY were known beforehand, one could use the following deterministic gradient updates to obtain θ * .EQUATION)Lemma 4. Let Assumption 2.3 be true. Then, choosingη k = O(k -a ) with 1/2 < a < 1, we have ∥ θ t -θ * ∥ = O exp(-t 1-a ) .Define the sequence δ t := θ t -θ t . We will establish the convergence rate of E ∥δ t ∥ 2 2 . From (13), and (27), we have the following expansion of δ t+1 .EQUATIONwhereQ t :=(I -α t+1 γ ⊤ * Σ Z γ * ), D t :=γ ⊤ * Σ Z γ * -γ t ⊤ Σ Z γ t .First we will establish an intermediate bound on E ∥δ t ∥ 2 . To do so, we will need the following result which shows that E ∥θ t -θ * ∥ 4 2 is bounded for all t which we prove in Section E.2. Lemma 5 (Boundedness of fourth moment of ∥θ t -θ * ∥). Let the conditions in Theorem 2 be true. Then, choosing α t , β t such that α t ≤ d -4κ-ϑ/2 z , and∞ t=1 (α 2 t + α t √ β t ) < ∞, we have E ∥θ t -θ * ∥ 4 2 is bounded by some constant M > 0. Lemma 6 (Intermediate bound on E[∥δ t ∥ 2 2 ]). Let the conditions in Theorem 2 be true. We have the following intermediate bound on E ∥δ t ∥ 2 :EQUATION)Proof of Lemma 6. Recall the update for δ t+1 obtained in (28).δ t+1 = Q t δ t + α t+1 D t θ t + α t+1 (γ t -γ * ) ⊤ Σ ZY -α t+1 γ t ⊤ ξ ZtYt + α t+1 γ t ⊤ ξ Zt γ t θ t .Then,EQUATIONThen, choosing α 1 (∥γ * ∥ 2 λ Z ) 2 < 1, using Young's inequality and Assumption 2.2, from (30) we get,E t ∥δ t+1 ∥ 2 2 ≤(1 -α t+1 µ)∥δ t ∥ 2 + 4α 2 t+1 ∥D t θ t ∥ 2 + ∥(γ t -γ * ) ⊤ Σ ZY ∥ 2 + 4α 2 t+1 ∥γ t ∥ 2 2 E ∥ξ ZtYt ∥ 2 2 + ∥γ t ∥ 4 2 E ∥ξ Zt ∥ 2 ∥θ t ∥ 2 + 2α t+1 δ t ⊤ Q t D t θ t + (γ t -γ * ) ⊤ Σ ZY ≲(1 -α t+1 µ)∥δ t ∥ 2 + 4α 2 t+1 ∥D t ∥ 2 ∥θ t ∥ 2 + ∥(γ t -γ * ) ⊤ Σ ZY ∥ 2 + 4Cα 2 t+1 d 2κ+ϑ/2 z + d 4κ+ϑ/2 z ∥θ t ∥ 2 + 2α t+1 δ t ⊤ Q t D t θ t + (γ t -γ * ) ⊤ Σ ZY ,where the last inequality follows by Assumption 3.1,and Assumption 3.3. Now, taking expectation on both sides, we obtainEQUATIONNow, the following bounds are true:1. We have thatEQUATIONwhere the first inequality follows by Cauchy-Schwarz inequality, the second inequality follows by (42), and Lemma 5.2. Using Σ ZY = O(1), and Lemma 3, we getEQUATION)3. We have thatEQUATIONwhere the first inequality follows by Hölder's inequality, the second inequality follows by Young's inequality, Cauchy-Schwarz inequality, and ∥Q t ∥ 2 < 1, and the third inequality follows by (42), and Lemma 5.Proof. Using the form of θ * , from (13) we get,EQUATIONwhereEQUATIONWe have the following bound on E ∥D t ∥ 4 2 by Lemma 3.E ∥D t ∥ 4 2 = E ∥(γ * -γ t ) ⊤ Σ Z γ * + γ t ⊤ Σ Z (γ * -γ t )∥ 4 2 = O(d 2+4κ z β 2 t ).(42) From (40), we haveEQUATIONEQUATIONEQUATIONThen, choosingEQUATION) We now have the following bounds:1. Using Assumption 3.1, and Assumption 3.3,EQUATION2. We have thatEQUATIONwhere, the first inequality follows by 3. We have that∥ Q t ∥ 2 = O(EQUATIONwhere, the first inequality follows by ∥ Q t ∥ 2 = O(EQUATION6. Using ∥ Q t ∥ 2 = O(1), Assumption 3.1, and Assumption 3.3,EQUATION)7. We have thatEQUATIONwhere, the first inequality follows by ∥ Q t ∥ 2 = O(1), and ∥Σ ZY ∥ 2 = O(1), the second inequality follows by Cauchy-Schwarz inequality, the third inequality follows by Lemma 3 and the fourth inequality follows by Young's inequality.8. Similar to (53), we have,EQUATION9. Using ∥ Q t ∥ 2 = O(1), Cauchy-Schwarz inequality, (50), and Young's inequality,EQUATION)10. By Assumption 2.2, we have,EQUATIONNow using Jensen's inequality, and combining (47), ( 48), ( 49), (50), and (51), we have,EQUATION)Combining ( 52), ( 53), (54), (55), and (56), we get,EQUATIONCombining ( 46), (57), and (58), we have,EQUATION) Now choosing α t , β t such that α t ≤ d -4κ-ϑ/2 z, andEQUATIONfor some constant 0 ≤ M < ∞.E.3 Comment on the convergence of (11)We now discuss the convergence properties of the update sequence (11), which we refer to as the conditional stochastic optimization (CSO) based updates, which we restate below:θ t+1 = θ t -α t+1 γ ⊤ t Z t (X ⊤ t θ t -Y t ), γ t+1 = γ t -β t+1 Z t (Z ⊤ t γ t -X ⊤ t ). Similar to (40), for the above updates, we have the following expansion:θ t+1 -θ * = Q t (θ t -θ * ) + α t+1 (γ t -γ * ) ⊤ Σ ZY + α t+1 D t θ * + α t+1 γ t ⊤ ξ Zt γ * (θ t -θ * ) + α t+1 γ t ⊤ ξ Zt γ * θ * + α t+1 γ t ⊤ ξ ZtYt -α t+1 γ t ⊤ Z t ϵ 2,t ⊤ θ t ,whereξ Zt = Σ Z -Z t Z ⊤ t , ξ ZtYt = Σ ZY -Z t Y t , Q t := I -α t+1 γ t ⊤ Σ Z γ * = Q t + α t+1 D t , and D t = (γ * -γ t )Σ Z γ * .Recall that the reason for the initial divergence of the updates in (11) are the potential negative eigenvalues of γ t ⊤ Σ Z γ * . Here we will show that if γ t ⊤ Σ Z γ * is positive semi-definite or γ t is close enough to γ * such that the negative eigenvalues (if any) are not too large in absolute values, then the updates in (11) indeed exhibit the same convergence rate as Algorithm 2. Assumption E.1. Let either of the following two conditions be true. For all t ≥ t 0 ,1. γ t Σ Z γ * is positive semidefinite. 2. ∥γ t -γ * ∥ 2 ≲ d z β t .Note that Condition 1 of Assumption E.1 is an idealized condition which is difficult to ensure for all t in reality. But of course if this is true, then γ t Σ Z γ * does not have a negative eigenvalue to cause divergence and the proof then follows exactly like Lemma 5.Hence, we will focus on the more realistic Condition 2 of Assumption E.1 which holds true almost surely [PJ92] . Since we are interested in the asymptotic rate of convergence of CSO updates (due to the requirement of Assumption E.1), we will only concentrate on the iterations t ≥ t 0 . In this case, the proof steps are similar to Theorem 2 except for two major differences, that we discuss below.Difference 1: Potential negative definiteness of γ t ⊤ Σ Z γ * :Under Condition 2, γ t ⊤ Σ Z γ * can indeed be negative definite. In general, if γ t ⊤ Σ Z γ * is negative definite then that is undesirable as we explain Section 3. In terms of the proof, we can no longer write(θ t -θ * ) ⊤ Q ⊤ t Q t (θ t -θ * ) ≤ ∥θ t -θ * ∥ 2 (which was possible to do in (43) in the proof of Lemma 5). Subsequently, (46) breaks down. But we will show that under Condition 2 the negative eigenvalues are not too large in terms of absolute values. Specifically, we can write,EQUATION)The term α 2 t+1 ∥D t ∥ 2 ∥θ t -θ * ∥ 2 is of the order of A 3,t defined in (45 t+1 ∥D t ∥) ≲ (1 + 2C γ α t+1 d κ z )which is too big. Recall that we want something at least of the order of α t+1 √ β t to show that θ t sequence is bounded. One could also try to use the fact that E [∥D t ∥] is small by Lemma 3. But since D t and θ t are interdependent, one needs to decouple them. One way to do this would be to use Cauchy-Shwarz inequalityas shown below.E ∥D t ∥∥θ t -θ * ∥ 2 ≤ E [∥D t ∥ 2 ] E [∥θ t -θ * ∥ 4 ] ≲ d z β t E [∥θ t -θ * ∥ 4 ].But that leads to the presence of E ∥θ t -θ * ∥ 4 in (43) which is potentially problematic due to the fact that on the left-hand side we have E ∥θ t+1 -θ * ∥ 2 .Difference 2: Presence of additional error term α t+1 γ t ⊤ Z t ϵ 2,t ⊤ θ t :When comparing (12) with ( 40), yet another crucial difference is the presence of the term α t+1 γ t ⊤ Z t ϵ 2,t ⊤ θ t . We will show by the following observations that this error term gets absorbed by other terms already present in (40) without affecting the convergence rate. Specifically, the following holds.1. Using the independence between Z, and ϵ 2,t , and by Assumption 2.2, we have,E t [( Q t (θ t -θ * ) + α t+1 (γ t -γ * ) ⊤ Σ ZY + α t+1 D t θ * + α t+1 γ t ⊤ ξ Zt γ * (θ t -θ * ) + α t+1 γ t ⊤ ξ Zt γ * θ * ) ⊤ γ t ⊤ Z t ϵ 2,t ⊤ θ t ] = 0.2. We also have thatα 2 t+1 E t (γ t ⊤ ξ ZtYt ) ⊤ γ t ⊤ Z t ϵ 2,t ⊤ θ t =α 2 t+1 (γ t ⊤ Σ Z γ t ∥θ * ∥ 2 + γ t ⊤ Σ Z γ t θ ⊤ * (θ t -θ * )) ≤α 2 t+1 (γ t ⊤ Σ Z γ t ∥θ * ∥ 2 + ∥γ t ⊤ Σ Z γ t (θ t -θ * )∥ 2 + ∥θ * ∥ 2 )This shows that the above term is of the same order as A 1,t and A 3,t defined in (44), and (45). 3. Finally, we haveα 2 t+1 E t ∥γ t ⊤ Z t ϵ 2,t ⊤ θ t ∥ 2 ≲ α 2 t+1 (∥γ t ∥ 2 ∥θ t -θ * ∥ 2 + ∥γ t ∥ 2 ∥θ * ∥ 2 ). So this term is of the order of A 3,t as well.Combining the above facts and following similar procedure as the proof of Theorem 2, one can show that the CSO updates achieve a similar rate under additional Assumption E.1.• It should be clear whether the error bar is the standard deviation or the standard error of the mean. Justification: Sections 2 and 3. Our work is mainly about theory and algorithm design, so there is no societal impact of the work performed. 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).Note that the streaming algorithm was not present in version 1, i.e.,[DVB23, v1]. 


Decision: Four reviewers have reviewed the paper, and their overall evaluation is positive. I concur with their assessment and believe the paper makes a strong contribution with compelling results. Two reviewers requested additional experiments and comparisons, to which the authors have responded satisfactorily. I strongly recommend that the authors include these additional simulations in the final version of the paper.

Reviews: review: paper_summary: paper_summary: This paper proposes and analyzes on-line algorithms for instrumental variable regression (IVaR) with streaming data.  Specifically, the authors consider the model:

$Y = g_{\theta^*}(X) + \epsilon_1$ 

where the covariate $X$ and noise $\epsilon_1$ are possibly correlated, but an instrumental variable $Z$ is available that satisfies:

$X = h_{\gamma^*}(Z) + \epsilon_2$ 

with $\epsilon_2$ being a centered (unobserved) noise.  Building on the prior works [MMLR20] and [HZCH20], the authors consider the IVaR problem formulated as a conditional stochastic optimization problem, presented in Equation (2):
 
$Minimize_g F(g) := E_Z E_{Y|Z} [ ( Y - E_{X|Z} [g(X)])^2 ].$

Considering a parameterized family of regression functions $G := \{ g(\theta; X) | \theta \}$, the authors observe that the gradient of F admits the expression as in Equation (3):

$\nabla F(\theta) = E_Z[ (E_{X|Z}[g(\theta;X)] - E_{Y|Z}[Y]) \cdot \nabla_{\theta} E_{X|Z}[g(\theta; X)]].$


This paper propose and analyze two streaming algorithms to solve the optimization problem formulation above. Notably, the proposals in this work do not require reformulating it to a minimax optimization problem as done in [MMLR20] or employing a nested sampling technique to reduce the bias in gradient estimation as in [HZCH20]. Specifically, the authors assume the availability of an oracle that can generate a sample $X$ (or two independent samples $X$ and $X’$) conditioned on $Z$, and then propose a stochastic gradient descent for IVaR. Additionally, the authors establish the rate of convergence assuming linear models. The proposed algorithms and claims are supported by numerical experiments with synthetic data.
strength_weakness: strength_weakness: While this work makes several significant theoretical contributions, especially in advancing the analysis, there are areas for potential improvement:

**1. Motivation for IVaR Problems (with Streaming Data)**: The importance of IVaR problems, especially with streaming data, should be highlighted more. Discussing example scenarios in addition to citing references would better motivate and convince readers of the problem's relevance.

**2. Further Experimental Validation**: Although this work is primarily theoretical, augmenting the experiment section with a more comprehensive set of experiments would be beneficial. It would be particularly valuable to see how the proposed algorithms perform on real-world datasets.
questions: questions: 1. While Algorithm 1 is generally applicable with the availability of two-sample oracles, Algorithm 2 seems to hinge critically on linear models. The authors note in lines 244-245 that a detailed treatment of the nonlinear case is left for future work due to its complexity in analysis, but I am curious if designing a working algorithm based on the insights in this paper would be feasible at least. Could the authors provide insights on the following: 

    (a) How would Algorithm 2 (or a variant) perform with nonlinear models?

    (b) Do the authors believe extending Algorithm 2 to nonlinear models is feasible, and what challenges might arise in the algorithm's design?

2. In line 57, can authors elaborate on why "Eq. (3) implies that one does not need the nested sampling technique to reduce the bias" with more details?

3. Suggestions

    (a) **Lines 150 - 151**:  I suggest the authors state Assumption 2.3 as follows: "... $P_{Y|X}$.  There exist constants such that $C_x, C_y, C_{xx}, C_{yx} > 0$ such that ..."

    (b) **Line 161**:  I guess it might make more sense to move Assumption 2.4 up to Line 147, right after the sentence "... pre-collected dataset."

    (c) **Page 9**: The authors may want to summarize the information in Lines 315 - 323 in the captions of Figures 3 and 4 for readers' convenience.  Also, it could be helpful to summarize the experimental results in the main text explicitly.

    (d) **Miscellaneous/Potential typos**: 

        i) Line 52: $X|Z$ instead of $Z|X$?

        ii) Line 154: remove "under"

        iii) Eq. (10) and Line 159: Please consider ending the sentence in Eq. (1) and start a new sentence in Line 159. Also, "if" in Line 159 should be capitalized.
limitations: limitations: This is primarily a theoretical work, and the authors discussed the potential limitations of their assumptions.
paper_summary: paper_summary: The paper shows algorithms for instrumental variable regression that dont need matirx inversions and mini-batches. At the same time, the paper give the rates of convergence.
strength_weakness: strength_weakness: None
questions: questions: I am curious about how the matrix inversion harm a algorithm performance in this setting.
limitations: limitations: None
paper_summary: paper_summary: This paper tackles the instrumental variable regression. The problem setting assumes a model $Y=g_{\theta^*}(X) + \epsilon_1$, but unlike the ordinary regression model, there are correlations between $X$ and $\epsilon_1$. The model assumes in addition an instrument variable $Z$ such that $Y$ and $X$ are independent conditional on $Z$, and $X=h_{\gamma^*}(Z)+\epsilon_2$. The target is to estimate $\theta^*$. The canonical approach is to use the two-stage least square (2SLS) method, where we first regress $X$ on $Z$ and then regress $Y$ on the $\widehat{X}$ that is estimated by $Z$. 

This paper considers an online setting where samples $(X_t, Y_t, Z_t)$ arrive sequentially. They propose two stochastic gradient descent (SGD) based algorithms. The first algorithm assumes that for each $Z_t$ we can resample $X_t'$ from the conditional distribution of $X | Z=Z_t$. The second algorithm does not make this assumption and can be considered an SGD version of the 2SLS algorithm. The paper provides the L2 convergence guarantee for both algorithms and they design experiments to demonstrate the effectiveness of their algorithms.
strength_weakness: strength_weakness: 1.  The contribution of the paper appears marginal. The first algorithm requires data resampling and this limits its applicability. The second algorithm is an intuitive adaption of the 2SLS algorithm to the SGD setting. The theorems are proved in the simple linear setting.
2. Apart from avoiding matrix inversion, the paper lacks necessary explanations as to why we should prefer an SGD version of 2SLS over the canonical offline 2SLS. It would be beneficial if more content were devoted to the advantages of Algorithm 2 compared to traditional benchmarks. For example, a theoretical comparison with 2SLS or an empirical comparison using real-world data would be helpful.
3. Could you explain why using two different synthetic settings for Algorithm 1 and Algorithm 2 in the Numerical Experiments section?
questions: questions: See "Weakness".
limitations: limitations: See "Weakness".
paper_summary: paper_summary: This paper proposes two stochastic optimization algorithms for instrumental variable regression (IVaR) that operate on streaming data without requiring matrix inversions or mini-batches. When the true model is linear, the paper proves that TOSG-IVaR converges at a rate of $\mathcal{O}(\log T/T)$ and OTSG-IVaR at a rate of $O(1/T^{1 - \iota})$ for any $\iota > 0$, where $T$ is the number of iterations. The proposed approaches avoid the "forbidden regression" problem of having to estimate the nuisance parameter relating $Z$ and $X$. Numerical experiments validate the theoretical results.
strength_weakness: strength_weakness: 1. The paper does not provide a comprehensive comparison with existing IVaR methods, particularly in terms of computational complexity and empirical performance on real-world datasets. Such a comparison would strengthen the paper's contributions. 

2. The assumptions made for the theoretical analysis, such as the identifiability conditions and moment assumptions, may be restrictive in some practical settings. A discussion on the robustness of the proposed algorithms under violations of these assumptions would be valuable.
questions: questions: The theoretical analysis focuses on the linear model setting. What are the challenges in extending the analysis to non-linear models, and what kind of assumptions or modifications would be required?
limitations: limitations: The limitations and future works are thoroughly discussed in the paper. Besides from that, I do not see more limitations.

Comments: title: Thank you for the insightful question
Response to the Authors' Rebuttal
comment: Dear Reviewer, 

Since the discussion period is about to end in a day, we wanted to make sure that our response fully addresses your question. If we have answered your question satisfactorily, we would appreciate if you consider increasing your score as you see fit.  

Thank you.
It is indeed true that learning  $\gamma_*$ requires incorporating $\nabla h$ in the update equation as we have alluded to in our global response in the section "**Alg. 2 convergence proof intuition**" but could not elaborate on due to lack of space. We elaborate on it here. 

**The same stability issue happens in this framework as we saw in the linear setting (Line 220 - 226) and is avoided by the same trick that we used in the linear setting.**
Consider the following version of the equation (11) adapted to the non-linear setting.
 $$\theta_{t+1}=\theta_t-\alpha_{t+1} h_{\gamma_t}(Z_t)(X_t^\top \theta_t-Y_t)\qquad
 \gamma_{t+1}=\gamma_t-\beta_{t+1} \nabla h_{\gamma_t}(Z_t)^\top( h_{\gamma_t}(Z_t)-X_t).\qquad \text{(11NL)}$$
Equation (11NL), similar to equation (12) in the main paper,  can be expanded in the following manner.

$$\theta_{t+1}-\theta_* = {\hat{Q}}\_t^{NL} ( \theta_t-\theta_* )+\alpha_{t+1}\mathbb{E}\_{\gamma_t}[(h_{\gamma_t}(Z_t)-h_{\gamma_*}(Z_t)) Y_t)]+\alpha_{t+1}D_t^{NL}\theta_*+\alpha_{t+1}\left(\mathbb{E}\_{\gamma_t}[h_{\gamma_t}(Z_t)h_{\gamma_*}(Z_t)^\top]-h_{\gamma_t}(Z_t)h_{\gamma_*}(Z_t)^\top\right)(\theta_t-\theta_*)$$ $$~~~~~~~~~~~~~~~~~~~~+\alpha_{t+1}\left(\mathbb{E}\_{\gamma_t}[h_{\gamma_t}(Z_t)h_{\gamma_*}(Z_t)^\top]-h_{\gamma_t}(Z_t)h_{\gamma_*}(Z_t)^\top\right)\theta_*+\alpha_{t+1}((h_{\gamma_t}(Z_t)-\mathbb{E}\_{\gamma_t}[h_{\gamma_t}(Z_t)])Y_t)
-\alpha_{t+1}h_{\gamma_t}(Z_t)\epsilon_{2,t}^\top\theta_t. \qquad \text{(12NL)}$$

where $\hat{Q}\_t^{NL}=(I-\alpha_{t+1}\mathbb{E}\_{\gamma_t}[h_{\gamma_t}(Z_t)h_{\gamma_*}(Z_t)^\top])$, $D_t^{NL}=\mathbb{E}\_{\gamma_t}[(h_{\gamma_*}(Z_t)-h_{\gamma_t}(Z_t)))h_{\gamma_*}(Z_t)^\top]$. 

First, let's focus on the **stability issue** associated with the first term on the RHS, i.e., $\hat{Q}\_t^{NL}(\theta_t-\theta_*)$. Just like the linear setting, here too, the matrix $\mathbb{E}\_{\gamma_t}[h_{\gamma_t}(Z_t)h_{\gamma_*}(Z_t)^\top]$ is not guaranteed to be positive semi-definite. So, we replace the term $X_t^\top\theta_t$ in equation (11NL) by $h_{\gamma_t}(Z_t)^\top\theta_t$ which leads to the following modified Algorithm 2 updates.
$$\theta_{t+1}=\theta_t-\alpha_{t+1} h_{\gamma_t}(Z_t)(h_{\gamma_t}(Z_t)^\top \theta_t-Y_t)\qquad \text{(13NL)}$$
$$\gamma_{t+1}=\gamma_t-\beta_{t+1} \nabla h_{\gamma_t}(Z_t)^\top( h_{\gamma_t}(Z_t)-X_t),\qquad \text{(14NL)}$$
where $\gamma_t\in\mathbb{R}^{d_\gamma}$, and $\nabla h_{\gamma_t}(Z_t)\in\mathbb{R}^{d_x\times d_\gamma}$ is the Jacobian of $h$ with respect to $\gamma_t$. 

For (13NL), we will have $\hat{Q}\_t^{NL}$ of the form $\hat{Q}\_t^{NL}=(I-\alpha_{t+1}\mathbb{E}\_{\gamma_t}[h_{\gamma_t}(Z_t)h_{\gamma_t}(Z_t)^\top])$. Here $\mathbb{E}\_{\gamma_t}[h_{\gamma_t}(Z_t)h_{\gamma_t}(Z_t)^\top]$ is positive semi-definite leading to the stability of the dynamics just like the linear case. 

**Now, we just have to show that rest of the terms on the right hand side of (12NL) converges similar to (12).** Recall that, except the first term, we control all the other terms on the right hand side of equation (12) mainly by using the martingale-difference property of the noise variables, and Lemma 3, i.e. the convergence of $\mathbb{E}[\||\gamma_t-\gamma_*\||^2]$. In (12NL), the martingale-difference property of the noise variables in the fourth to seventh terms on the right hand side clearly holds.

It remains to obtain a result analogous to Lemma 3. To do so, we look at equation (14NL). Analysis of equation (14NL) to establish a convergence rate for the $\gamma_t$ updates to $\gamma_*$ or $h_{\gamma_t}$ to $h_{\gamma_*}$ is straightforward as long as $H(\gamma)\coloneqq\mathbb{E}[\||X-h_{\gamma}(Z)\||^2]$ is strongly-convex [PJ92], or satisfies Polyak-Łojasiewicz (PL) inequality [KNS16]. This increases the model flexibility considerably as PL inequalities are satisfied for a wide class of non-linear DNN models [LZB22]. Beyond strong-convex and PL cases, the analysis is challenging although the same algorithmic framework still applies as a methodology.

Putting the above pieces together, it is possible to obtain the rates of convergence for the case when $Z$ to $X$ is non-linear.

Please reach out to us if you have any additional question.
I thank the authors for addressing my concerns and questions.  I find most of their responses satisfactory, and would like to encourage them to incorporate the explanations and remarks into the revision.  

However, I am still not clearly seeing how Algorithm 2 can be immediately applied to the scenario where the relationship between $X$ and $Z$ is non-linear, as the authors claim.  Specifically, the update rules in Eqs. (13) and (14) -- which stem from Eq. (11) and include an additional modification trick to promote stability as discussed in Lines 220 -- 226 -- seem to rely on the linear model assumption.  Thus, it is unclear to me how these would translate to the non-linear model setting described by the authors, namely, the setting where $Y = {\theta^*}^{\top} X + \epsilon_1$ with $X = h_{\gamma^*}(X) + \epsilon_2$.  I believe the update rule should involve $\nabla h_{\gamma_t}$ and possibly a similar modification to promote algorithmic stability.

Could you please clarify the update rules for this setting by specifying the counterparts of Eqs. (13) and (14), and any necessary modifications (if needed) that corresponds to replacing $g(\theta_t, X_t) = X_t^{\top} \theta_t$ with $Z_t^{\top} \gamma_t \theta_t$ in the linear setting?
Thank you for the detailed reply, I am now more convinced of the paper's technical contribution. I will raise my score from 5 to 6.
Thank you very much for your detailed response. I will keep my original score.