Differentially Private Two-Stage Gradient Descent for Instrumental Variable Regression

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.

50,401 characters · 11 sections · 43 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.

abstractWe study instrumental variable regression (IVaR) under differential privacy constraints. Classical IVaR methods (like two-stage least squares regression) rely on solving moment equations that directly use sensitive covariates and instruments, creating significant risks of privacy leakage and posing challenges in designing algorithms that are both statistically efficient and differentially private. We propose a noisy two-stage gradient descent algorithm that ensures $\rho$-zero-concentrated differential privacy by injecting carefully calibrated noise into the gradient updates. Our analysis establishes finite-sample convergence rates for the proposed method, showing that the algorithm achieves consistency while preserving privacy. In particular, we derive precise bounds quantifying the trade-off among optimization, privacy, and sampling error. To the best of our knowledge, this is the first work to provide both privacy guarantees and provable convergence rates for instrumental variable regression in linear models. We further validate our theoretical findings with experiments on both synthetic and real datasets, demonstrating that our method offers practical accuracy-privacy trade-offs.

Introduction

Instrumental variable regression (IVaR) is a foundational tool in causal inference, designed to recover structural parameters when standard estimators fail due to endogeneity. In many observational settings, covariates are influenced by unobserved confounders, causing naive methods (such as the ordinary least squares (OLS) in the context of linear regression) to produce biased and inconsistent estimates. IVaR circumvents this by leveraging instruments, which are variables that are predictive of the endogenous regressors but independent of hidden confounders, to enable consistent estimation of causal effects hausman2001mismeasured,wooldridge2010econometric,10.1257/jep.15.4.69. This perspective is increasingly important in machine learning, for example in recommendation systems where user exposure is confounded by prior preferences si2022model, or in reinforcement learning where actions and rewards are jointly influenced by unobserved context xu2023instrumental. In such settings, IVaR provides a principled way to disentangle causal effects from spurious correlations, enabling more reliable decision making.

However, many applications of IVaR involve sensitive data, such as individual health records, financial transactions, or user interactions, where protecting privacy is of paramount importance. In such settings, releasing model estimates or even intermediate statistics can leak information about individuals in the dataset. Differential privacy (DP) dwork2006calibrating provides a mathematically rigorous framework to ensure that an algorithm’s output does not reveal sensitive information about any single data point. Despite the importance of IVaR in causal inference, to the best of our knowledge, there are no prior works addressing the problem of performing IVaR under differential privacy. This gap motivates the central question of this paper:

center[center omitted — 159 chars of source]

Our work focuses on answering this question in the context of linear IVaR models. To situate our contributions, we briefly review existing work on DP methods for OLS regression, with additional discussion in Section (ref). Several predominant approaches have emerged in the literature: (i) perturbion methods, where the empirical covariance and cross-covariance matrices are privatized before solving the normal equations; (ii) consensus-based methods, including propose-test-release and exponential mechanism approaches, which directly privatize the estimator through carefully designed randomized output rules; and (iii) gradient perturbation methods, where iterative optimization algorithms are made private by clipping gradients and injecting calibrated Gaussian noise. While all three approaches ensure differential privacy, gradient perturbation combined with clipping has been shown to yield the sharpest statistical rates in OLS regression, particularly in high-dimensional and finite-sample regimes bassily2014private, brown2024private.

Given the centrality of IVaR in causal inference, it is natural to explore whether the aforementioned techniques can be adapted to this setting. Unlike OLS, however, IVaR is based on moment conditions involving both covariates and instruments, making it less straightforward to design private algorithms. In particular, sufficient-statistics perturbation and consensus-based methods have not been explored, and their adaptation is non-trivial due to the inherent ill-posedness of IVaR under weak instruments and the sensitivity of the moment equations. These challenges stem from the fact that the closed-form 2SLS estimator depends on nested matrix multiplications and inversions, whose sensitivities are difficult to characterize directly. To address this, we design a two-stage gradient-descent-based algorithm that enables injecting noise at the iteration level in a principled manner. This structure allows us to rigorously control sensitivity, calibrate the noise, and derive non-asymptotic utility guarantees.

Specifically, we make the following contributions in this work:

itemize[noitemsep, leftmargin=-0.001in] • We introduce DP-2S-GD (Algorithm (ref)), the first differentially private algorithm for instrumental variable regression, built on noisy gradient descent with gradient clipping. • We establish finite-sample non-asymptotic error bound for DP-2S-GD (Theorem (ref)), explicitly characterizing the trade-off between optimization, privacy, and sampling error. The main technical challenge is to carefully control the interaction between privacy-induced noise and the contraction of the gradient dynamics across iterations, with the privacy guarantee ensured by Proposition (ref). • We validate our theoretical analysis with experiments on synthetic and real-world datasets, demonstrating practical accuracy-privacy trade-offs (Section (ref)).

By developing differentially private algorithm for IVaR and establishing its theoretical bound, we enable researchers to perform robust causal analyses without compromising the confidentiality of individuals in the dataset.

Related work

Differential Privacy for Regression. One can group private regression methods into the following broad families. (1) Output/objective perturbation (private empirical risk minimization (ERM)): add noise to the final estimator (output perturbation) or inject a random linear/quadratic term into a strongly convex loss before optimizing (objective perturbation); these one-shot mechanisms give $(\varepsilon, \delta)$-DP guarantees and excess-risk bounds for convex ERM (chaudhuri2011differentially; kifer2012private; bassily2014private). Recent refinements, e.g. redberg2023improving, leverage subsampling and tighter accounting to improve accuracy. (2) Sufficient-statistics (matrix) perturbation: release noisy surrogates of ($\mathbf{X}^{\top}\mathbf{X}, \mathbf{X}^{\top}\mathbf{y}$) (or related second-moment structures) and then solve the (regularized) normal equations; this route enables OLS-specific inference but can suffer under ill-conditioning because noise is injected at the Gram-matrix level (dwork2014analyze; sheffet2017differentially). tsfadia2022friendlycore proposes a subsample-and-aggregate framework that can, in principle, be adapted to regression settings by privately estimating the relevant sufficient statistics on carefully selected data subsets. Further developments in this direction include bernstein2019differentially and ferrando2024private. (3) Exponential mechanism: privately selects an output by randomly choosing among candidates with probabilities that grow exponentially with their quality score, with parameters controlling how strongly it favors the higher-scoring options. This mechanism is frequently applied in constructing algorithm to privately select a regression model from a pool of non-private OLS fits on subsets of the data (ramsay2021differentially, cumings2022differentially, amin2022easy). (4) Gradient perturbation (DP-(S)GD): clip per-example (mini-batch or full) gradients and add Gaussian noise at each step, tracking privacy with bounded log moment generating function of privacy loss random variable wang2019subsampled, R\'enyi DP, and subsampled-RDP-which yields tight composition for many small releases and scales well to large $n, p$ without forming $\mathbf{X}^{\top} \mathbf{X}$. (abadi2016deep; bun2016concentrated; mironov2017renyi; wang2019subsampled).

We favor gradient perturbation for multi-stage estimators like IVaR because it (i) composes tightly across many noisy steps using modern privacy accountants, (ii) avoids spectrum-dependent blow-ups from noising $\mathbf{X}^{\top} \mathbf{X}$ (sheffet2017differentially) and (iii) yields strong convergence rates while fitting standard training pipelines (including using minibatches, streaming, early stopping) and enabling modular, stage-wise design, which is preferable for practice (bassily2014private, abadi2016deep). Although there exists DP techniques for estimating gram matrices that avoid blow-ups, e.g., via carefully calibrated noise or regularization (brown2023fast, kamath2019privately), the purely sufficient statistics-based pipelines require larger sample sizes (polynomial to condition number) than gradient-based approaches to reach a comparable accuracy in high dimensions {brown2024insufficient}. That said, we note that the convergence of private first-order gradient methods still depend heavily on the condition number varshney2022nearly, liu2023near.

Instrumental Variable Regression (IVaR) has been extensively studied in econometrics 10.1257/jep.15.4.69,angrist2009mostly. Classical methods such as two-stage least squares (2SLS) admit closed-form solutions but face limitations in modern applications: they do not scale well to high-dimensional or streaming data, cannot easily incorporate regularization, and are restricted to linear models. This has motivated optimization-based approaches, including convex–concave formulations of nonlinear IV muandet2020dualinstrumentalvariableregression, stochastic optimization methods for scalable and online estimation della2023stochastic,chen2024stochasticoptimizationalgorithmsinstrumental,peixoto2024nonparametric, and bi-level gradient descent algorithms with convergence guarantees liang2025transformers. Extensions to nonlinear IV include kernel-based methods singh2019kernel and DeepIV hartford2017deep. Despite these advances, prior work assumes unrestricted access to the data and does not provide end-to-end differential privacy guarantees, which are increasingly critical in sensitive domains such as healthcare, finance, and online platforms. To our knowledge, no existing method offers DP guarantees with finite-sample convergence rates for linear IV/2SLS that explicitly account for instrument strength, sample size, dimension, and iteration complexity.

Notations:\quad Throughout this paper, unless otherwise specified, we use lower-case letters to denote random variable or individual data samples, and upper-case letters to denote datasets, i.e. collections of samples. Bolded letters represent vectors and matrices, whereas unbolded letters represent scalars.

Preliminaries

Privacy notions

We first review widely used notions of privacy in the literature. Two datasets $D$ and $D'$ are said to be neighbors if they differ in exactly one entry. The concept of neighboring datasets allows us to formally quantify the level of differential privacy. The two most common notions are $(\varepsilon,\delta)$-differential privacy and zero-concentrated differential privacy (zCDP).

defn[$(\varepsilon,\delta)$-Differential Privacy dwork2006calibrating] A randomized mechanism $M$ satisfies $(\varepsilon,\delta)$-differential privacy if for all neighboring datasets $D,D'$ and all measurable sets $S$, we have $\Pr[M(D) \in S] \leq e^{\varepsilon} \Pr[M(D') \in S] + \delta.$ Here $\varepsilon \geq 0$ controls the multiplicative privacy loss, while $\delta \in [0,1]$ allows for a small probability of arbitrary deviation.

defn[Zero-Concentrated Differential Privacy (zCDP) dwork2016concentrated,bun2016concentrated] A randomized mechanism $M$ satisfies $\rho$-zero-concentrated differential privacy ($\rho$-zCDP) if for all neighboring datasets $D,D'$ and all $\alpha > 1$, we have the $D_{\alpha}\!\left(M(D)\,\|\,M(D')\right) \le \rho \alpha$, where $D_{\alpha}(P\|Q)$ denotes the R\'enyi divergence (see Appendix (ref) for the definition) of order $\alpha$ between distributions $P$ and $Q$.

While $(\varepsilon,\delta)$-DP is the most widely used notion of privacy, it can be too coarse for analyzing iterative mechanisms, as composition accumulates $\varepsilon$ and $\delta$ linearly\footnote{$(\epsilon,\delta)$-DP also satisfies an advanced composition that asymptotically matches zCDP, but the composition is more cumbersome and typically less practical than the clean additive composition offered by zCDP.}. In contrast, zero-concentrated differential privacy (zCDP) characterizes privacy loss through R\'enyi divergences, which ensures that the privacy loss random variable enjoys a sub-Gaussian concentration property. This yields two key benefits: (i) tighter composition, since zCDP parameters add under composition, and (ii) smooth conversion, since $\rho$-zCDP implies $(\varepsilon,\delta)$-DP with $\varepsilon = \rho + 2\sqrt{\rho \log(1/\delta)}$; see bun2016concentrated. As a result, we choose zCDP for technical convenience since it provides simple additive composition rule and leaner formulas in our context where we compose a large number of identical Gaussian mechanisms across both stages in 2SLS algorithm.

IVaR Model and Assumptions

Endogeneity is a central challenge in linear regression. Suppose we aim to estimate the causal effect of the regressor $\mathbf{x}\in\mathbb{R}^p$ on the outcome $y\in\mathbb{R}$. However, there exists an unobserved confounder $\mathbf{u}$ that affects both $\mathbf{x}$ and $y$, thereby violating the standard exogeneity assumption that $\mathbf{x}$ is uncorrelated with the noise. As a result, the OLS estimator becomes biased and inconsistent. Instrumental variable regression (IVaR) is a widely adopted method to handle endogeneity by including $\mathbf{z}\in\mathbb{R}^q$, an instrumental variable (IV), to the model 10.1257/jep.15.4.69:

equation[equation omitted — 223 chars of source]

where the error terms $\epsilon_1$ and $\boldsymbol{\epsilon}_2$ are correlated due to the common confounder $\boldsymbol{u}$; see Figure (ref) for an illustration. Given the dataset $(\mathbf{Z}, \mathbf{X}, \mathbf{Y})=\left\{\left(\mathbf{z}_i, \mathbf{x}_i, y_i\right)\right\}_{i=1}^{n}$\footnote{Throughout this paper, we assume each entry of the dataset is independently and identically distributed (i.i.d.).}, the objective of the IVaR model is to solve the following bi-level optimization problem:

equation[equation omitted — 477 chars of source]

Optimization problem (ref) admits a closed-form solution. A classical approach to solve (ref) is the two-stage least squares (2SLS) estimator; see Definition (ref).

defn[2SLS estimator] Given observational data $(\mathbf{Z},\mathbf{X},\mathbf{Y})=\{(\mathbf{z}_i,\mathbf{x}_i, y_i)\}_{i=1}^n$, the 2SLS estimator $\hat{\boldsymbol{\beta}}_{\textsf{2SLS}}$ is obtained through two consecutive OLS regressions: \begin{itemize}[noitemsep] • First stage: Regress $\mathbf{X}$ on $\mathbf{Z}$ to obtain $\hat{\boldsymbol{\Theta}}$ \begin{align*} \hat{\boldsymbol{\Theta}} = (\mathbf{Z}^\top\mathbf{Z})^{-1}\mathbf{Z}^\top\mathbf{X}. \end{align*} • Second stage: Regress $\mathbf{Y}$ on $\hat{\mathbf{X}}:=\mathbf{Z}\hat{\boldsymbol{\Theta}}$ to obtain $\hat{\boldsymbol{\beta}}_{\textsf{2SLS}}$: \begin{align*} \hat{\boldsymbol{\beta}}_{2SLS} = (\hat{\boldsymbol{\Theta}}^\top\mathbf{Z}^\top\mathbf{Z}\hat{\boldsymbol{\Theta}})^{-1}\hat{\boldsymbol{\Theta}}^\top\mathbf{Z}^\top\mathbf{Y}. \end{align*} \end{itemize}

In the following sections, we will use $\hat{\boldsymbol{\beta}}$ to denote the 2SLS estimator for simplicity. We impose the following standard assumptions for IVaR model.

asp[IVaR Assumptions] A random variable $\mathbf{z}\in \mathbb{R}^q$ is a valid IV, if it satisfies: \begin{itemize}[noitemsep] • Fully identification: $q\geq p$ (without loss of generality, we assume data $\mathbf{Z},\mathbf{X}$ are full rank). • Correlation to $\boldsymbol{x}$: $\textsf{Corr}(\mathbf{z},\mathbf{x})\neq\mathbf{0}$. • Exclusion to $y$: $\textsf{Corr}(\mathbf{z},\epsilon_1)=\mathbf{0}$. \end{itemize}

In Assumption (ref), condition (i) ensures the existence of the unique solution $\hat{\boldsymbol{\beta}}$ in (ref), condition (ii) guarantees that the instrument explains nontrivial variation in the endogenous regressor $\mathbf{x}$, and condition (iii) ensures that the instrument affects the outcome $y$ only through $\mathbf{x}$. These conditions are crucial for eliminating endogeneity and achieving consistent estimation for $\boldsymbol{\beta}$. See stock2011introduction for a detailed discussion.

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

We further impose the following assumptions to establish non-asymptotic rates.

aspWe assume the following conditions hold: \begin{itemize}[noitemsep] • $\mathbf{z}$ is a mean-zero isotropic sub-Gaussian random vector. That is, $\mathbb{E}[\mathbf{z}]=\mathbf{0}$, $\mathbb{E}[\mathbf{z}\mathbf{z}^\top]=\mathbf{I}_q$, and for some $\sigma_z>0$, $\mathbb{E}[e^{{u\langle\mathbf{z}_{i},\mathbf{v}\rangle}}]\leq \exp\{\frac{u^2\sigma_z^2\|\mathbf{v}\|^2}{2}\},\forall u\in\mathbb{R}, \mathbf{v}\in\mathbb{R}^q$. • $\epsilon_1,\boldsymbol{\epsilon}_2$ are mean-zero sub-Gaussian. That is, $\mathbb{E}[\epsilon_1]=0, \mathbb{E}[\boldsymbol{\epsilon}_2]=\mathbf{0}$, and for some $\sigma_1,\sigma_2>0$, $\mathbb{E}[e^{u \epsilon_1}]\leq \exp\{\frac{u^2\sigma_1^2}{2}\}$, and $\mathbb{E}[e^{u \langle\boldsymbol{\epsilon}_2,\mathbf{v}\rangle}]\leq \exp\{\frac{u^2\sigma_2^2\|\mathbf{v}\|^2}{2}\},\forall u\in\mathbb{R},\mathbf{v}\in\mathbb{R}^p$. \end{itemize}

Assumption (ref) provides the minimal conditions required to leverage concentration results from high-dimensional random design analysis Vershynin_2018. Specifically, with condition (i), we have the high-probability concentration bound for the empirical covariance matrix $\frac{\mathbf{Z}^\top\mathbf{Z}}{n}$ (see Lemma (ref)). Condition (ii) further ensures high-probability concentration of the cross terms $\frac{\mathbf{Z}^\top\boldsymbol{\mathcal{E}}_1}{n}$ and $\frac{\mathbf{Z}^\top\boldsymbol{\mathcal{E}}_2}{n}$ (see Lemma (ref)), where $(\boldsymbol{\mathcal{E}}_1,\boldsymbol{\mathcal{E}}_2)=\{(\epsilon_{1,i},\boldsymbol{\epsilon}_{2,i})\}_{i=1}^n$ denotes the sample realization of errors. With these conditions, we derive high-probability concentration bound for the sample covariance matrix of $\hat{\mathbf{X}}:=\mathbf{Z}\hat{\boldsymbol{\Theta}}$ (see Lemma (ref)), and finally establish the non-asymptotic error bound $\|\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}\|$ (see Lemma (ref)).

Privacy in IVaR may be required at different levels depending on the application. In some cases, protecting only the causal effect $\boldsymbol{\beta}$ is sufficient, for instance when the first-stage compliance relation $\boldsymbol{\Theta}$ is public, secondary, or not sensitive. In other cases, privacy must also extend to the first-stage parameter $\boldsymbol{\Theta}$, such as when instruments involve sensitive behavioral data, proprietary mechanisms, or institutional policies. To ensure end-to-end privacy in the IVaR model, we adopt the framework of zCDP. We allocate two privacy parameters: $\rho_1$ for the first-stage parameter estimates $\{\boldsymbol{\Theta}^{(t)}\}_{t=1}^T$, and $\rho_2$ for the second-stage parameter estimates $\{\boldsymbol{\beta}^{(t)}\}_{t=1}^T$. By the composition property of zCDP, the overall procedure satisfies $(\rho_1+\rho_2)$-zCDP.

Algorithm and Theoretical Guarantees

We begin with a baseline two-stage gradient descent algorithm, denoted as 2S-GD, for solving the IVaR problem (ref). The detailed procedure is deferred to Appendix (ref), Algorithm (ref). The method alternates between two coupled updates at each iteration: (i) updating the first-stage projection matrix $\boldsymbol{\Theta}^{(t)}$, which maps instruments $\mathbf{Z}$ to covariates $\mathbf{X}$, and (ii) updating the second-stage regression parameter $\boldsymbol{\beta}^{(t)}$ based on the predicted covariates. This iterative procedure can be viewed as a gradient-based analogue of the classical two-stage least squares estimator.

In this section, we propose a differentially private two-stage gradient descent algorithm, termed DP-2S-GD, to solve the IVaR problem (ref) while ensuring rigorous privacy guarantees. The algorithm is summarized in Algorithm (ref). Compared with 2S-GD, DP-2S-GD incorporates two key modifications: (i) per-sample clipping is applied to gradients in both stages to bound the sensitivity of each update, ensuring that no single datapoint can disproportionately affect the results, and (ii) Gaussian perturbations are injected into both the $\boldsymbol{\Theta}$- and $\boldsymbol{\beta}$-updates at every iteration, with noise scales calibrated to the target privacy budgets $\rho_1$ and $\rho_2$.

algorithm[algorithm omitted — 1,518 chars of source]

The privacy analysis proceeds by treating the two stages as separate Gaussian mechanisms with sensitivity controlled by clipping parameters $\gamma_1$ and $\gamma_2$. By the properties of zero-concentrated differential privacy, the choice of noise scales $\lambda_1,\lambda_2$ uniquely determines the effective privacy losses $\rho_1,\rho_2$, which compose additively across iterations. Consequently, for any pre-specified privacy budgets $(\rho_1,\rho_2)$, one can calibrate $(\lambda_1,\lambda_2)$ to ensure that DP-2S-GD achieves the desired privacy guarantees. We next establish formal theoretical results, including both privacy accounting and utility bounds for the resulting estimators.

propIf we set $\lambda_1 = \frac{2\gamma_{1}}{n}\sqrt{\frac{T}{\rho_1}}$ and $\lambda_2 = \frac{2\gamma_{2}}{n}\sqrt{\frac{T}{\rho_2}}$, Algorithm (ref) is $\rho$-zCDP, where $\rho:=\rho_1+\rho_2=\frac{2 T}{n^2}\left(\frac{\gamma_1^2}{\lambda_1^2}+\frac{\gamma_2^2}{\lambda_2^2}\right)$.

The proof of Proposition (ref) is provided in Appendix (ref).

rmkProposition (ref) highlights several tradeoffs among the parameters. To preserve the same privacy levels $\rho_1, \rho_2$, the noise scales $\lambda_1, \lambda_2$ must increase with larger clipping thresholds $\gamma_1, \gamma_2$, or with larger number of iterations $T$. Conversely, a larger sample size $n$ allows for smaller noise scales while maintaining the same privacy guarantees.

thmFor any fixed $\boldsymbol{\Theta}\in\mathbb{R}^{q\times p}$ and $\boldsymbol{\beta}\in\mathbb{R}^{p}$, consider the Algorithm (ref) with fixed step sizes satisfying \begin{align} 0<\eta<\frac{2}{(1+\delta(\tau))^2}, \quad 0<\alpha<\frac{4}{2\bar{\gamma}(\tau)+\gamma(\tau)}, \end{align} under Assumption (ref), with parameters \begin{equation} \begin{aligned} \lambda_1 = \frac{2\gamma_{1}}{n}\sqrt{\frac{T}{\rho_1}},\quad\lambda_2 = \frac{2\gamma_{2}}{n}\sqrt{\frac{T}{\rho_2}},\quad \gamma_1 = \gamma_2 = c_0\left(\sqrt{q}+\sqrt{\tau+\log(nT)}\right)^2, \end{aligned} \end{equation} and number of iterations \begin{align} T\lesssim \frac{\rho_1 n^{2-\epsilon}}{p(\sqrt{q}+\sqrt{\tau})^6}, \end{align} where $\epsilon>0$ is a small constant. If \begin{align} n\geq c_1\max\left\{pq(\tau+\log(pq))^2, \frac{\left(\sqrt{q}+\sqrt{\tau}\right)^3}{\sqrt{\min\{\rho_1,\rho_2\}}}\right\}, \end{align} for any fixed $\tau$, with probability $1-c_2e^{-\tau}$, we have \begin{align} \begin{split} \|\boldsymbol{\beta}^{(T)}-\hat{\boldsymbol{\beta}}\| &\lesssim \kappa(\tau)^{\frac{T}{2}}+\frac{\sqrt{p}(\sqrt{q}+\sqrt{\tau})^3}{n\sqrt{\min\{\rho_1,\rho_2\}}}\sqrt{T}+\frac{\sqrt{pq}(\tau+\log(pq))}{\sqrt{n}}, \end{split} \end{align} where $0<\kappa(\tau)<1$ is the contraction rate, $\delta(\tau)>0$ is a numerically small term, and $\bar{\gamma}(\tau), \underline{\gamma}(\tau)$ are the high-probability upper/lower bounds on the eigenvalues of $\frac{\hat{\boldsymbol{\Theta}}^\top\mathbf{Z}^\top\mathbf{Z}\hat{\boldsymbol{\Theta}}}{n}$. The specific definitions of $\delta(\tau), \bar{\gamma}(\tau), \underline{\gamma}(\tau),$ and $\kappa(\tau)$ are deferred to (ref).

The proof of Theorem (ref) is presented in Appendix (ref). We now offer several remarks regarding this theorem. In the presentation of Theorem (ref), all constants $c_0, c_1, c_2$ and scaling factors hidden in "$\lesssim$" are independent of major parameters $n, p, q, T, \rho_1, \rho_2, \tau$. These constants only depend on problem-specific parameters $\boldsymbol{\beta}, \boldsymbol{\Theta}, \sigma_z, \sigma_1, \sigma_2$.

rmkConsider the population optimization problem $ \min_{\boldsymbol{\beta}}\tilde{\mathcal{L}}(\boldsymbol{\beta})=\mathbb{E}\left[(y-\mathbf{z}^\top\boldsymbol{\Theta\beta})^2\right]$, and the (deterministic) two-stage gradient descent algorithm: \begin{align*} \boldsymbol{\Theta}^{(t+1)}=\boldsymbol{\Theta}^{(t)}-\eta_{GD}\mathbb{E}\left[\mathbf{z}(\mathbf{z}^\top\boldsymbol{\Theta}^{(t)}-\mathbf{x}^\top)\right],\quad\quad \boldsymbol{\beta}^{(t+1)}=\boldsymbol{\beta}^{(t)}-\alpha_{GD}\mathbb{E}\left[\boldsymbol{\Theta}^{\top}\mathbf{z}(\mathbf{z}^\top\boldsymbol{\Theta}\boldsymbol{\beta}^{(t)}-y)\right]. \end{align*} It can be easily shown that under Assumption (ref), the sufficient condition for learning rates to guarantee monotonic convergence are \begin{align*} 0<\eta_{GD}<2,\quad 0<\alpha_{GD} < \frac{2}{\|\boldsymbol{\Theta}\|^2}. \end{align*} We note that in our learning rate condition (ref), we introduce $\delta(\tau)$ and $\psi(\tau)$ to account for the randomness in data. If we have infinite samples, the condition (ref) becomes \begin{align*} 0<\eta<2,\quad 0<\alpha<\frac{4}{2\|\boldsymbol{\Theta}\|^2+\sigma_{\min}^2(\boldsymbol{\Theta})}. \end{align*} Comparing to $\eta_{GD}$ and $\alpha_{GD}$, notice that we have the same $\eta$ condition. However, the $\alpha$ condition is slightly tighter to control the randomness introduced by the first-stage estimates $\boldsymbol{\Theta}^{(t)}$.

rmkFrom (ref), the optimal contraction rate $\kappa^\star(\tau)$ is achieved when the learning rates are set as \begin{align} \eta_{approx}^\star = \frac{2}{(1+\delta(\tau))^2+(1-\delta(\tau))^2},\quad \alpha_{approx}^\star = \frac{2}{\bar{\gamma}(\tau)+\gamma(\tau)}. \end{align} In this case, we have \begin{gather*} \kappa_{\boldsymbol{\beta}}^\star(\tau) = \frac{\bar{\gamma}(\tau)}{\bar{\gamma}(\tau)+\gamma(\tau)},\quad \kappa_{\boldsymbol{\Theta}}^\star(\tau) = \frac{(1+\delta(\tau))^2-(1-\delta(\tau))^2}{(1+\delta(\tau))^2+(1-\delta(\tau))^2}, \kappa^\star(\tau) = \max\left\{\kappa_{\boldsymbol{\beta}}^\star(\tau), \kappa_{\boldsymbol{\Theta}}^\star(\tau)\right\}. \end{gather*} We emphasize that although $\eta_{\text{approx}}^\star$ and $\alpha_{\text{approx}}^\star$ minimize the contraction rate, they should be viewed as approximately optimal step sizes, as the scaling constants in the bound (ref) vary with different choices of step sizes. Empirically, we also observe that the estimator's error is fairly insensitive in a neighborhood of $\eta_{\text{approx}}^\star, \alpha_{\text{approx}}^\star$. See Appendix (ref) for empirical results.

rmkFrom Proposition (ref), the choice of $\lambda_1, \lambda_2$ in (ref) guarantees that Algorithm (ref) is $\rho$-zCDP. The parameters $\gamma_1$ and $\gamma_2$ are selected so that, with high probability, the clipping operation does not alter the gradients; see Lemma (ref) for details.

rmkThe error bound (ref) consists of three dominant terms. The first term $\kappa(\tau)^{\frac{T}{2}}$ characterizes the convergence of the gradient descent algorithm, which decays exponentially with $T$. The second term $\frac{\sqrt{p}(\sqrt{q}+\sqrt{\tau})^3}{n\sqrt{\min\{\rho_1,\rho_2\}}}\sqrt{T}$ captures the cumulative effect of the injected Gaussian noise, which grows with $\sqrt{T}$ due to the parameter choices in (ref) that ensure privacy. The third term $\frac{\sqrt{pq}(\tau+\log(pq))}{\sqrt{n}}$ represents the inherent statistical error in estimating $\hat{\boldsymbol{\beta}}$ via noiseless gradient descent, which decreases with larger sample size $n$. This decomposition highlights the trade-offs between convergence phase and privacy requirement, while also accounting for the structural statistical accuracy attainable from gradient descent.

rmkThe condition for $T$ in (ref) is necessary to control the noise scale $\lambda_1$ in Proposition (ref), since the derivation of (ref) relies on the high-probability concentration of $\|\boldsymbol{\Theta}^{(T)}-\hat{\boldsymbol{\Theta}}\|$. With limited sample size $n$, if $\rho_1$ is small, i.e. we want high privacy on $\boldsymbol{\Theta}^{(1)}, \ldots, \boldsymbol{\Theta}^{(T)}$, we can only set a moderate number of iterations $T$, otherwise the bound (ref) doesn't hold. See Section (ref) for experiments.

rmkFor given sample size $n$, the dominating terms for each $T$ range are: \begin{equation*} \|\boldsymbol{\beta}^{(T)}-\hat{\boldsymbol{\beta}}\| \;\lesssim\; \begin{cases} \kappa(\tau)^{\tfrac{T}{2}}, & if T \leq \dfrac{\log\!\left(\tfrac{n}{pq(\tau+\log(pq))^2}\right)} {\log\!\left(\tfrac{1}{\kappa(\tau)}\right)}, \\[1.2em] \dfrac{\sqrt{pq}\,(\tau+\log(pq))}{\sqrt{n}}, & if \dfrac{\log\!\left(\tfrac{n}{pq(\tau+\log(pq))^2}\right)} {\log\!\left(\tfrac{1}{\kappa(\tau)}\right)} < T \leq \dfrac{n \min\{\rho_1,\rho_2\}q(\tau+\log(pq))^2}{(\sqrt{q}+\sqrt{\tau})^6}, \\[1.2em] \dfrac{\sqrt{p}(\sqrt{q}+\sqrt{\tau})^3}{n\sqrt{\min\{\rho_1,\rho_2\}}}\; \sqrt{T}, & if \dfrac{n \min\{\rho_1,\rho_2\}q(\tau+\log(pq))^2}{(\sqrt{q}+\sqrt{\tau})^6} < T \lesssim \dfrac{\rho_1 n^{2-\epsilon}}{p(\sqrt{q}+\sqrt{\tau})^6}. \end{cases} \end{equation*}

Hence, the optimum number of iterations $T$ is sub-linear but super-logarithmic to $n$. Figure (ref) qualitatively illustrates the trend of the error bound (ref) as a function of $T$. This is consistent with our experimental observations in Section (ref).

figure[figure omitted — 219 chars of source]

coroConsider running Algorithm (ref) with $\rho_1=\infty$ and $\rho_2=\infty$ (i.e. no privacy provided). For any $T>0$, the bound (ref) is dominated by \begin{align} \|\boldsymbol{\beta}^{(T)}-\hat{\boldsymbol{\beta}}\| &\lesssim \kappa(\tau)^{\frac{T}{2}}+\frac{\sqrt{pq}(\tau+\log(pq))}{\sqrt{n}}, \end{align} which is exactly the convergence rate of the $\texttt{2S-GD}$ algorithm (ref).

rmkWe note that when $\rho_1,\rho_2=\infty$ and $T\rightarrow\infty$, the error rate (ref) still has an additional $\sqrt{p}$ factor compared to the error rate of 2SLS estimator $\|\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}\|$ (see Lemma (ref) for the precise statement). This performance gap is an inherent limitation of gradient-based approximations to 2SLS, which is further confirmed by simulations in Appendix (ref). As an intuitive explanation, the closed-form 2SLS estimator solves both stages in (ref) using sample moments, whereas the gradient-descent procedure has to approximate the second-stage moment condition. In the population, the optimality condition for 2SLS ensures $\mathbb{E}[\mathbf{z}(\mathbf{y}-\mathbf{x}^\top\boldsymbol{\beta}^\star)]=0$. However, in the finite-sample gradient-descent iteration, the update direction will eventually depend on $\tfrac{1}{n}\mathbf{Z}^\top(\mathbf{Y}-\mathbf{Z}\hat{\boldsymbol{\Theta}}\hat{\boldsymbol{\beta}}),$ which involves the empirical residual $\mathbf{r}:=\mathbf{Y}-\mathbf{Z}\hat{\boldsymbol{\Theta}}\hat{\boldsymbol{\beta}}$. In particular, the quantity $\|\tfrac{1}{n}\mathbf{Z}^\top\mathbf{r}\|$ has an error rate $\tfrac{\sqrt{pq}(\tau+\log(pq))}{\sqrt{n}}$ (see Lemma (ref)), which yields the additional $\sqrt{p}$ factor.

rmkIn practice, the intermediate estimates $\{\boldsymbol{\Theta}^{(t)}\}_{t=1}^T$ are not always required to be released, so in some settings it suffices to ensure privacy only for $\{\boldsymbol{\beta}^{(t)}\}_{t=1}^T$. In Algorithm (ref), setting $\rho_1=\infty$ implies that no noise $\boldsymbol{\Xi}^{(t)}$ needs to be injected in the first stage, and we can simply return $\{\boldsymbol{\beta}^{(t)}\}_{t=1}^T$ under privacy budget $\rho_2$. Under this regime, the error bound (ref) continues to hold, except that the condition on $T$ in (ref) is no longer required. See Appendix (ref) for further details.

Experiments

We conduct experiments using both synthetic data and real data to validate our theoretical findings. For all experiments, we set $\tau=5$, and step sizes $\eta=\frac{1}{(1+\delta(\tau))^2}$, $\alpha=\frac{2}{2\bar{\gamma}(\tau)+\underline{\gamma}(\tau)}$. As a practical guideline, $\rho=0.1$ is considered as strong privacy, $\rho=1$ is considered as moderate privacy, and $\rho=10$ is considered as weak privacy\footnote{The corresponding $(\epsilon,\delta)$-DP values using the conversion formula $\epsilon = \rho + 2\sqrt{\rho \log(1/\delta)}$ (with $\delta = 10^{-5}$): $\rho = 0.1 \Leftrightarrow (\epsilon,\delta) = (2.25,10^{-5})$, $\rho = 1 \Leftrightarrow (\epsilon,\delta) = (7.79,10^{-5})$, and $\rho = 10 \Leftrightarrow (\epsilon,\delta) = (31.47,10^{-5})$.}.

Synthetic Data Simulations

We generate synthetic data according to the IVaR model in (ref). To simulate the correlation between $\boldsymbol{\epsilon}_1$ and $\epsilon_2$, we include a confounder $\mathbf{u}\in\mathbb{R}^r$, and set $\epsilon_1 = \boldsymbol{\Phi}^\top\mathbf{u}_i+\boldsymbol{\epsilon}_{x}$ and $\epsilon_2 = \boldsymbol{\phi}^\top\mathbf{u}+\epsilon_{y}$, and generate each entry of the dataset $(\mathbf{Z}, \mathbf{X}, \mathbf{Y})=\{(\mathbf{z}_i,\mathbf{x}_i,y_i)\}_{i=1}^n$ according to the following model: $\mathbf{x}_i=\boldsymbol{\Theta}^\top \mathbf{z}_i+\boldsymbol{\Phi}^\top\mathbf{u}_i+\boldsymbol{\epsilon}_{x,i},$ and $\mathbf{y}_i=\boldsymbol{\beta}^\top\mathbf{x}_i+\boldsymbol{\phi}^\top\mathbf{u}_i+\epsilon_{y,i},$ where the ground-truth parameters are $\boldsymbol{\beta}\in\mathbb{R}^p, \boldsymbol{\Theta}\in\mathbb{R}^{q\times p}$, $\mathbf{\Phi}\in\mathbb{R}^{r\times p}$, $\boldsymbol{\phi}\in\mathbb{R}^r$. These parameters are drawn as follows: $\boldsymbol{\beta}\sim \mathcal{N}(\mathbf{0},\mathbf{I}_p)$, $\boldsymbol{\Theta}\sim 5\mathbf{I}_{q\times p}+\mathbf{E}$ with $\mathbf{E}_{ij}\sim\mathcal{N}(0,1)$. $\boldsymbol{\Phi}_{ij}\sim \mathcal{N}(0,1)$, and $\boldsymbol{\phi}\sim\mathcal{N}(\mathbf{0},\mathbf{I}_r)$. For each simulation, we then sample $\mathbf{z}_i\sim \mathcal{N}(\mathbf{0},\mathbf{I}_q)$, $\mathbf{u}_i\sim \mathcal{N}(\mathbf{0},\mathbf{I}_r)$, $\boldsymbol{\epsilon}_{x,i}\sim \mathcal{N}(\mathbf{0},\mathbf{I}_p)$, and $\epsilon_{y,i}\sim \mathcal{N}(0,1)$.

figure[figure omitted — 859 chars of source]

Figure (ref) compares the performance of Algorithm (ref) across different sample sizes $n$ under varying privacy allocations. We fix the total privacy budget at $\rho=\rho_1+\rho_2=10$, set the number of iterations to $T=20$, and examine three regimes: (i) $\rho_1=1, \rho_2=9$, (ii) $\rho_1=5, \rho_2=5$, and (iii) $\rho_1=9, \rho_2=1$. In Figure (ref)(a), with $p=q=r=5$, all points lie in the plateau region of Figure (ref), so the error decreases at the rate $\tfrac{1}{\sqrt{n}}$. In contrast, Figure (ref)(b) sets $p=q=r=50$. Here, $T=20$ violates condition (ref), leading to significantly larger errors compared to Figure (ref)(a).

figure[figure omitted — 758 chars of source]

The impact of $T$ is further investigated in Figure (ref), from which we observe that, with limited sample size $n$, if we enforce high privacy guarantee on $\{\boldsymbol{\Theta}^{(t)}\}_{t=1}^T$ (i.e. with small $\rho_1$), the error grows significantly after certain $T$ is reached. This cutoff aligns with the condition on $T$ specified in (ref). In contrast, when privacy is required only for $\{\boldsymbol{\beta}^{(t)}\}_{t=1}^T$ (i.e., with small $\rho_2$), the error behavior closely matches the theoretical predictions illustrated in Figure (ref).

Real-Data Experiments

We further evaluate our algorithm on the Angrist dataset laborsupply, which has been widely applied in the IVaR literature. This study examines the causal effect of children bearing on female labor supply, leveraging the gender composition of the first two children as an instrument\footnote{Research shows that parents whose first two children are of the same sex are significantly more likely to have an additional child Westoff1972. At the same time, the sex composition of the first two children can be treated as randomly assigned and is not directly related to the mother’s labor supply.}. The endogenous regressor $\mathbf{x}$ is the number of children bearing, the outcome $\mathbf{y}$ is the mother's labor supply measured in number of working weeks per year, and the instrument $\mathbf{z}$ is a binary variable indicating whether the first two children are of the same gender. The original dataset contains $394,835$ samples. For illustration purpose, we randomly draw a subset of $20,000$ samples and keep $n=8065$ effective observations with number of children $\geq 2$. We center all variables $\mathbf{z}, \mathbf{x}, \mathbf{y}$ and run Algorithm (ref) with $T=20$ iterations. Figure (ref) presents the results over 1000 independent runs with privacy budgets $\rho_1=1, \rho_2=1$. As shown in Figure (ref), the estimated $\boldsymbol{\beta}^{(T)}$ concentrates around $-4.3$, indicating that having an additional child reduces the mother's labor supply by approximately 4.3 weeks per year. This estimate is consistent with the 2SLS benchmark.

figure[figure omitted — 908 chars of source]

From Figure (ref), we observe that Algorithm (ref) converges in expectation after about 15 iterations. The dispersion of the estimates is determined by the privacy budgets: increasing $\rho_1$ and $\rho_2$ yield estimates that are more tightly concentrated around the 2SLS benchmark, while smaller budgets result in greater variability. Additional experiments are provided in Appendix (ref).

Conclusion

We have introduced DP-2S-GD, a differentially private two-stage gradient descent method for IVaR problem. The algorithm achieves $(\rho_1+\rho_2)$-zCDP by injecting carefully calibrated Gaussian noise. We have established finite-sample convergence guarantees that capture the trade-offs among optimization dynamics, privacy constraints, and statistical error. Our theoretical analysis shows that setting the number of iterations $T$ to be sub-linear yet super-logarithmic in $n$ minimizes the estimation error, a result that is corroborated by our experiments. We have further illustrated the practical utility of our method through an application to the Angrist dataset. On the other hand, we note that, regardless of the privacy constraint, the convergence of the two-stage gradient descent estimator to $\hat{\boldsymbol{\beta}}$ is slower by a $\sqrt{p}$ compared to the convergence of $\hat{\boldsymbol{\beta}}$ to the true parameter $\boldsymbol{\beta}$ (see Remark (ref)). Improving this rate (via algorithmic modifications) and establishing lower-bounds for privacy-accuracy tradeoffs for the IVaR problem are interesting future directions.

Acknowledgments

Krishnakumar Balasubramanian was supported in part by NSF grant DMS-2413426. Haodong Liang and Lifeng Lai were supported in part by NSF grants CCF-2232907, ECCS-2514514 and ECCS-2448268.