The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
50,401 characters
Differentially Private Two-Stage Gradient Descent for Instrumental Variable Regression
\maketitle
\begin{abstract}
We study \emph{instrumental variable regression} (IVaR) under \emph{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 \emph{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.
\end{abstract}
\section{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 \emph{instruments}, which are variables that are predictive of the endogenous regressors but independent of hidden confounders, to enable consistent estimation of causal effects \citep{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~\citep{si2022model}, or in reinforcement learning where actions and rewards are jointly influenced by unobserved context~\citep{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) \citep{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 \emph{no prior works} addressing the problem of performing IVaR under differential privacy.
This gap motivates the central question of this paper:
\begin{center}
\emph{Can we design differentially private algorithms for instrumental variable models \\that achieve statistically efficient convergence rates?}
\end{center}
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{sec:relwork}.
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 \citep{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 \textbf{contributions} in this work:\vspace{-0.1in}
\begin{itemize}[noitemsep, leftmargin=-0.001in]
\item We introduce \texttt{DP-2S-GD} (Algorithm \ref{alg: DP-2S-GD-II}), the first differentially private algorithm for instrumental variable regression, built on noisy gradient descent with gradient clipping.
\item We establish finite-sample non-asymptotic error bound for \texttt{DP-2S-GD} (Theorem \ref{thm: main result II}), 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{lem: privacy II - main}.
\item We validate our theoretical analysis with experiments on synthetic and real-world datasets, demonstrating practical accuracy-privacy trade-offs (Section \ref{sec: experiments}).
\end{itemize}\vspace{-0.1in}
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.
\subsection{Related work}\label{sec:relwork}
\textbf{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 (\cite{chaudhuri2011differentially}; \cite{kifer2012private}; \cite{bassily2014private}). Recent refinements, e.g. \cite{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 (\cite{dwork2014analyze}; \cite{sheffet2017differentially}). \cite{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~\cite{bernstein2019differentially} and ~\cite{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 (\cite{ramsay2021differentially}, \cite{cumings2022differentially}, \cite{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 \cite{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}$. (\cite{abadi2016deep}; \cite{bun2016concentrated}; \cite{mironov2017renyi}; \cite{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}$ (\cite{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 (\cite{bassily2014private}, \cite{abadi2016deep}). Although there exists DP techniques for estimating gram matrices that avoid blow-ups, e.g., via carefully calibrated noise or regularization (\cite{brown2023fast}, \cite{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 \citep{{brown2024insufficient}}. That said, we note that the convergence of private first-order gradient methods still depend heavily on the condition number \cite{varshney2022nearly, liu2023near}.
\textbf{Instrumental Variable Regression (IVaR)} has been extensively studied in econometrics~\citep{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~\citep{muandet2020dualinstrumentalvariableregression}, stochastic optimization methods for scalable and online estimation~\citep{della2023stochastic,chen2024stochasticoptimizationalgorithmsinstrumental,peixoto2024nonparametric}, and bi-level gradient descent algorithms with convergence guarantees~\citep{liang2025transformers}. Extensions to nonlinear IV include kernel-based methods~\citep{singh2019kernel} and DeepIV~\citep{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.
\textit{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.
\section{Preliminaries}
\subsection{Privacy notions}
We first review widely used notions of privacy in the literature. Two datasets $D$ and $D'$ are said to be \emph{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).
\begin{defn}[$(\varepsilon,\delta)$-Differential Privacy \citep{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.
\end{defn}
\begin{defn}[Zero-Concentrated Differential Privacy (zCDP) \citep{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{sec:renyidef} for the definition) of order $\alpha$ between distributions $P$ and $Q$.
\end{defn}
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) \emph{tighter composition}, since zCDP parameters add under composition, and (ii) \emph{smooth conversion}, since $\rho$-zCDP implies $(\varepsilon,\delta)$-DP with $\varepsilon = \rho + 2\sqrt{\rho \log(1/\delta)}$; see~\citet[Proposition 1.3]{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.
\subsection{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 \citep{10.1257/jep.15.4.69}:
\begin{equation}\label{eq: ivar model}
\begin{aligned}
y=\boldsymbol{\beta}^{\top}\mathbf{x}+\epsilon_1,\quad \quad
\mathbf{x}=\boldsymbol{\Theta}^\top\mathbf{z}+\boldsymbol{\epsilon}_2,
\end{aligned}
\end{equation}
where the error terms $\epsilon_1$ and $\boldsymbol{\epsilon}_2$ are correlated due to the common confounder $\boldsymbol{u}$; see Figure~\ref{fig: ivarfigure} 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:
\begin{equation}\label{eq: 2SLS problem}
\hat{\boldsymbol{\beta}}=\arg\min_{\boldsymbol{\beta}\in \mathbb{R}^p}\Big\{\mathcal{L}(\boldsymbol{\beta})=\frac{1}{n} \sum_{i=1}^n\left(y_i-\boldsymbol{\beta}^\top\hat{\boldsymbol{\Theta}}^\top\mathbf{z}_i \right)^2\Big\},
\text{s.t. }\hat{\boldsymbol{\Theta}}=\underset{\boldsymbol{\Theta}\in\mathbb{R}^{q\times p} }{\arg \min }\Big\{\frac{1}{n} \sum_{j=1}^n\|\mathbf{x}_j- \boldsymbol{\Theta}^\top \mathbf{z}_j\|^2\Big\}.
\end{equation}
Optimization problem \ref{eq: 2SLS problem} admits a closed-form solution. A classical approach to solve \eqref{eq: 2SLS problem} is the two-stage least squares (2SLS) estimator; see Definition~\ref{def:2sls}.
\begin{defn}[\textsf{2SLS} estimator] \label{def:2sls}
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:\vspace{-2mm}
\begin{itemize}[noitemsep]
\item [i.] \emph{First stage}: Regress $\mathbf{X}$ on $\mathbf{Z}$ to obtain $\hat{\boldsymbol{\Theta}}$ \vspace{-1mm}
\begin{align*}
\hat{\boldsymbol{\Theta}} = (\mathbf{Z}^\top\mathbf{Z})^{-1}\mathbf{Z}^\top\mathbf{X}.
\end{align*}\vspace{-4mm}
\item [ii.] \emph{Second stage}: Regress $\mathbf{Y}$ on $\hat{\mathbf{X}}:=\mathbf{Z}\hat{\boldsymbol{\Theta}}$ to obtain $\hat{\boldsymbol{\beta}}_{\textsf{2SLS}}$:\vspace{-1mm}
\begin{align*}
\hat{\boldsymbol{\beta}}_{\textsf{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}
\end{defn}\vspace{-2mm}
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.
\begin{asp}[IVaR Assumptions]\label{asp: IV}
A random variable $\mathbf{z}\in \mathbb{R}^q$ is a valid IV, if it satisfies:
\begin{itemize}[noitemsep]
\item[(i)] Fully identification: $q\geq p$ (without loss of generality, we assume data $\mathbf{Z},\mathbf{X}$ are full rank).
\item [(ii)] Correlation to $\boldsymbol{x}$: $\textsf{Corr}(\mathbf{z},\mathbf{x})\neq\mathbf{0}$.
\item [(iii)] Exclusion to $y$: $\textsf{Corr}(\mathbf{z},\epsilon_1)=\mathbf{0}$.
\end{itemize}
\end{asp}
In Assumption \ref{asp: IV}, condition (i) ensures the existence of the unique solution $\hat{\boldsymbol{\beta}}$ in~\eqref{eq: 2SLS problem}, 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 \citet[Chapter 12]{stock2011introduction} for a detailed discussion.
\begin{figure}[t]
\centering
\begin{tikzpicture}[
node distance=18mm and 23mm,
every node/.style={font=\small},
latent/.style={ellipse, draw, minimum width=26mm, minimum height=9mm},
observed/.style={rectangle, draw, minimum width=28mm, minimum height=9mm, rounded corners=2pt},
confound/.style={ellipse, draw, dashed, minimum width=28mm, minimum height=9mm, fill=gray!15},
>=Stealth
]
\node[latent] (Z) {Instrument $\mathbf{z}$};
\node[observed, right=of Z] (X) {Endogenous regressor $\mathbf{x}$};
\node[observed, right=of X] (Y) {Outcome $y$};
\node[confound, below=of $(X)!0.5!(Y)$] (U) {Unobserved confounder $\mathbf{u}$};
\draw[->, blue, thick] (Z) -- (X) node[midway, above] {\small relevance};
\draw[->, blue, thick] (X) -- (Y) node[midway, above] {\small causal effect};
\draw[->, dashed] (U) -- (X);
\draw[->, dashed] (U) -- (Y);
\draw[->, red, thick, dashed]
($(Z.east)+(0.01cm,0)$)
.. controls +(-1cm,1cm) and +(-1.4cm,1.8cm) ..
($(Y.west)+(0,0)$)
node[midway, allow upside down] (crosspos) {};
\node[font=\footnotesize\bfseries, red] at (crosspos) {X};
\node[font=\itshape, red, above=7mm of Y, xshift=-60mm] (note)
{Exclusion: No direct $\mathbf{z}\!\to\!y$ path};
\end{tikzpicture}
\caption{IVaR model: Instrument $\mathbf{z}$ is correlated with the endogenous regressor $\mathbf{x}$ and influences the outcome $y$ only indirectly through $\mathbf{x}$, while an unobserved confounder $\mathbf{u}$ affects both $\mathbf{x}$ and $y$.}
\label{fig: ivarfigure}
\end{figure}
We further impose the following assumptions to establish non-asymptotic rates.
\begin{asp}\label{asp: Z}
We assume the following conditions hold:
\begin{itemize}[noitemsep]
\item[(i)] $\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$.
\item[(ii)] $\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}
\end{asp}
Assumption \ref{asp: Z} provides the minimal conditions required to leverage concentration results from high-dimensional random design analysis \citep{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{lem: Z bound}). 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{lem: ZE2 bound}), 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{lem: lambda_min H star}), and finally establish the non-asymptotic error bound
$\|\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}\|$ (see Lemma \ref{lem: beta hat error bound}).
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.
\section{Algorithm and Theoretical Guarantees}\label{sec: Theoretical Guarantees}
We begin with a baseline two-stage gradient descent algorithm, denoted as \texttt{2S-GD}, for solving the IVaR problem \eqref{eq: 2SLS problem}. The detailed procedure is deferred to Appendix \ref{sec:renyidef}, Algorithm~\ref{alg: 2S-GD}. 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 \texttt{DP-2S-GD}, to solve the IVaR problem \eqref{eq: 2SLS problem} while ensuring rigorous privacy guarantees. The algorithm is summarized in Algorithm~\ref{alg: DP-2S-GD-II}. Compared with \texttt{2S-GD}, \texttt{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$.
\begin{algorithm}[t]
\caption{DP-2S-GD}\label{alg: DP-2S-GD-II}
\begin{algorithmic}[1]
\State \textbf{Input:} Data $\mathbf{Z}\in\mathbb{R}^{n\times q}$, $\mathbf{X}\in\mathbb{R}^{n\times p}$, $\mathbf{Y}\in\mathbb{R}^n$, target privacy budgets $\rho_1, \rho_2>0$, step sizes $\eta, \alpha > 0$, number of iterations $T$
\State \textbf{Parameters:} Noise scales $\lambda_1, \lambda_2>0$, clipping thresholds $\gamma_1, \gamma_2 > 0$
\State Initialize $\boldsymbol{\beta}^{(0)} =\mathbf{0}_p$, $\boldsymbol{\Theta}^{(0)} =\mathbf{0}_{q \times p}$
\For{$t = 0,1,\ldots,T-1$}
\State Draw $\boldsymbol{\Xi}^{(t)}$ with $\mathrm{vec}(\boldsymbol{\Xi}^{(t)})\sim \mathcal{N}(\mathbf{0},\lambda_{1}^{2}\mathbf{I}_{q}\otimes\mathbf{I}_{p})$
\State Draw $\boldsymbol{\nu}^{(t)}\sim \mathcal{N}(\mathbf{0}, \lambda_{2}^2\mathbf{I}_{p})$
\State $\boldsymbol{\Theta}^{(t+1)}=\boldsymbol{\Theta}^{(t)}-\frac{\eta}{n}\sum_{i=1}^n\text{CLIP}_{\gamma_1}\!\left(\mathbf{z}_i(\mathbf{z}_i^\top\boldsymbol{\Theta}^{(t)}-\mathbf{x}_i^\top)\right)+\eta\boldsymbol{\Xi}^{(t)}$
\State $\boldsymbol{\beta}^{(t+1)}=\boldsymbol{\beta}^{(t)}-\frac{\alpha}{n}\sum_{i=1}^{n}\text{CLIP}_{\gamma_2}\!\left(\boldsymbol{\Theta}^{(t)\top}\mathbf{z}_i(\mathbf{z}_{i}^{\top}\boldsymbol{\Theta}^{(t)}\boldsymbol{\beta}^{(t)}-y_{i})\right)+\alpha\boldsymbol{\nu}^{(t)}$
\EndFor
\State \Return $\{\boldsymbol{\Theta}^{(t)}\}_{t=1}^T,\ \{\boldsymbol{\beta}^{(t)}\}_{t=1}^T$
\end{algorithmic}
\end{algorithm}
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 \texttt{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.
\begin{prop}\label{lem: privacy II - main}
If 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{alg: DP-2S-GD-II} 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)$.
\end{prop}
The proof of Proposition \ref{lem: privacy II - main} is provided in Appendix \ref{sec: proof of lem: privacy II}.
\begin{rmk}
Proposition \ref{lem: privacy II - main} 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.
\end{rmk}
\begin{thm}\label{thm: main result II}
For any fixed $\boldsymbol{\Theta}\in\mathbb{R}^{q\times p}$ and $\boldsymbol{\beta}\in\mathbb{R}^{p}$, consider the Algorithm \ref{alg: DP-2S-GD-II} with fixed step sizes satisfying
\begin{align}\label{eq: learning rates condition}
0<\eta<\frac{2}{(1+\delta(\tau))^2}, \quad 0<\alpha<\frac{4}{2\bar{\gamma}(\tau)+\underline{\gamma}(\tau)},
\end{align}
under Assumption \ref{asp: Z}, with parameters
\begin{equation}\label{eq: parameter settings}
\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}\label{eq: T condition thm 1}
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}\label{eq: n condition}
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}\label{eq: e_beta(T) bound thm2}
\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 \eqref{eq: definitions of parameters}.
\end{thm}
The proof of Theorem \ref{thm: main result II} is presented in Appendix \ref{sec: proof of main result II}.
We now offer several remarks regarding this theorem. In the presentation of Theorem \ref{thm: main result II}, 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$.
\begin{rmk}\label{rmk: comparison to gd}
Consider 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{asp: Z}, the sufficient condition for learning rates to guarantee \emph{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 \eqref{eq: learning rates condition}, we introduce $\delta(\tau)$ and $\psi(\tau)$ to account for the randomness in data. If we have infinite samples, the condition \eqref{eq: learning rates condition} 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)}$.
\end{rmk}
\begin{rmk}
From \eqref{eq: definitions of parameters}, the optimal contraction rate $\kappa^\star(\tau)$ is achieved when the learning rates are set as
\begin{align}\label{eq: optimum step sizes}
\eta_{\text{approx}}^\star = \frac{2}{(1+\delta(\tau))^2+(1-\delta(\tau))^2},\quad \alpha_{\text{approx}}^\star = \frac{2}{\bar{\gamma}(\tau)+\underline{\gamma}(\tau)}.
\end{align}
In this case, we have
\begin{gather*}
\kappa_{\boldsymbol{\beta}}^\star(\tau) = \frac{\bar{\gamma}(\tau)}{\bar{\gamma}(\tau)+\underline{\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 \eqref{eq: e_beta(T) bound thm2} 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{sec: exp step size tuning} for empirical results.
\end{rmk}
\begin{rmk}
From Proposition \ref{lem: privacy II - main}, the choice of $\lambda_1, \lambda_2$ in \eqref{eq: parameter settings} guarantees that Algorithm \ref{alg: DP-2S-GD-II} 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{lem: no clipping II} for details.
\end{rmk}
\begin{rmk}
The error bound \eqref{eq: e_beta(T) bound thm2} 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 \eqref{eq: parameter settings} 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.
\end{rmk}
\begin{rmk}
The condition for $T$ in \eqref{eq: parameter settings} is necessary to control the noise scale $\lambda_1$ in Proposition~\ref{lem: privacy II - main}, since the derivation of \eqref{eq: e_beta(T) bound thm2} 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 \eqref{eq: e_beta(T) bound thm2} doesn't hold. See Section \ref{sec: experiments} for experiments.
\end{rmk}
\begin{rmk}
For 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}},
& \text{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}},
& \text{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},
& \text{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*}
\end{rmk}
Hence, the optimum number of iterations $T$ is sub-linear but super-logarithmic to $n$. Figure \ref{fig: sketch} qualitatively illustrates the trend of the error bound \eqref{eq: e_beta(T) bound thm2} as a function of $T$. This is consistent with our experimental observations in Section \ref{sec: experiments}.
\begin{figure}[t]
\centering
\includegraphics[width=0.6\textwidth]{Figures/sketch.png}
\caption{Qualitative trend of the error bound \eqref{eq: e_beta(T) bound thm2} as a function of $T$.}
\label{fig: sketch}
\end{figure}
\begin{coro}\label{coro: main result}
Consider running Algorithm \ref{alg: DP-2S-GD-II} with $\rho_1=\infty$ and $\rho_2=\infty$ (i.e. no privacy provided). For any $T>0$, the bound \eqref{eq: e_beta(T) bound thm2} is dominated by
\begin{align}\label{eq: e_beta(T) bound ordinary}
\|\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{alg: 2S-GD}.
\end{coro}
\begin{rmk}\label{rmk:compto2sls}
We note that when $\rho_1,\rho_2=\infty$ and $T\rightarrow\infty$, the error rate \eqref{eq: e_beta(T) bound ordinary} still has an additional $\sqrt{p}$ factor compared to the error rate of 2SLS estimator $\|\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}\|$ (see Lemma \ref{lem: beta hat error bound} 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{sec: exp convergence rate compare}. As an intuitive explanation, the closed-form 2SLS estimator solves both stages in \eqref{eq: 2SLS problem} 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{lem: Zr bound}), which yields the additional $\sqrt{p}$ factor.
\end{rmk}
\begin{rmk}
In 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{alg: DP-2S-GD-II}, 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 \eqref{eq: e_beta(T) bound thm2} continues to hold, except that the condition on $T$ in \eqref{eq: T condition thm 1} is no longer required. See Appendix \ref{sec: privacy for beta only} for further details.
\end{rmk}
\section{Experiments}\label{sec: 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})$.}.
\subsection{Synthetic Data Simulations}\label{sec: synthetic experiments}
We generate synthetic data according to the IVaR model in \eqref{eq: ivar model}. 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)$.
\begin{figure}[t]
\centering
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics[width=\linewidth]{Figures/Alg1_vs_n_p=5_q=5.png}
\caption{}
\end{subfigure}
\hfill
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics[width=\linewidth]{Figures/Alg1_vs_n_p=50_q=50.png}
\caption{}
\end{subfigure}
\caption{Comparison of Algorithm \ref{alg: DP-2S-GD-II}'s performance versus $n$. We set $T=20$, (a) $p=q=5$, (b) $p=q=50$. Note that the $T$ condition \eqref{eq: parameter settings} is not satisfied in (b). We set the total budget $\rho=10$ and compare three regimes: (i) $\rho_1=1, \rho_2=9$, (ii) $\rho_1=5, \rho_2=5$, (iii) $\rho_1=9, \rho_2=1$. The curves are averaged over 100 runs, with vertical bars representing the standard errors.}
\label{fig: error_vs_n}
\end{figure}
Figure \ref{fig: error_vs_n} compares the performance of Algorithm \ref{alg: DP-2S-GD-II} 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{fig: error_vs_n}(a), with $p=q=r=5$, all points lie in the plateau region of Figure \ref{fig: sketch}, so the error decreases at the rate $\tfrac{1}{\sqrt{n}}$. In contrast, Figure \ref{fig: error_vs_n}(b) sets $p=q=r=50$. Here, $T=20$ violates condition \eqref{eq: T condition thm 1}, leading to significantly larger errors compared to Figure \ref{fig: error_vs_n}(a).
\begin{figure}[t]
\centering
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics[width=\linewidth]{Figures/Alg1_vs_T_n=1000_p=5_q=5_rho1_small.png}
\caption{}
\end{subfigure}
\hfill
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics[width=\linewidth]{Figures/Alg1_vs_T_n=1000_p=5_q=5_rho2_small.png}
\caption{}
\end{subfigure}
\caption{Comparison of Algorithm \ref{alg: DP-2S-GD-II}'s performance versus number of iterations $T$. We fix $n=1000$, $p=q=5$, (a) keep $\rho_2$ large and vary $\rho_1$, (b) keep $\rho_1$ large and vary $\rho_2$. The curves are averaged over 100 runs, with vertical bars representing the standard errors.}
\label{fig: error_vs_T}
\end{figure}
The impact of $T$ is further investigated in Figure \ref{fig: error_vs_T}, 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 \eqref{eq: T condition thm 1}. 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{fig: sketch}.
\subsection{Real-Data Experiments}
We further evaluate our algorithm on the Angrist dataset \citep{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 \citep{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{alg: DP-2S-GD-II} with $T=20$ iterations. Figure \ref{fig:Angrist_results_rho1=1_rho2=1} presents the results over 1000 independent runs with privacy budgets $\rho_1=1, \rho_2=1$. As shown in Figure \ref{fig:Angrist_boxplot_rho1=1_rho2=1}, 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.
\begin{figure}[t]
\centering
\begin{subfigure}{0.48\linewidth}
\centering
\includegraphics[width=0.95\textwidth]{Figures/Angrist_Boxplot_rho1=1_rho2=1_T=20.png}
\subcaption{}
\label{fig:Angrist_boxplot_rho1=1_rho2=1}
\end{subfigure}
\hfill
\begin{subfigure}{0.48\linewidth}
\centering
\includegraphics[width=0.9\textwidth]{Figures/Angrist_Parameter_Path_rho1=1_rho2=1_T=20.png}
\subcaption{}
\label{fig:Angrist_parameter_path_rho1=1_rho2=1}
\end{subfigure}
\caption{Results on the Angrist dataset with $T=20, \rho_1=1, \rho_2=1$.
(a) Boxplot of estimated $\boldsymbol{\beta}^{(T)}$, over 1000 runs. (b) Learning paths of parameters $\boldsymbol{\beta}^{(t)}, \boldsymbol{\Theta}^{(t)}$, over 1000 runs. The shaded area represents the standard error.\vspace{-0.21in}}
\label{fig:Angrist_results_rho1=1_rho2=1}
\end{figure}
From Figure \ref{fig:Angrist_parameter_path_rho1=1_rho2=1}, we observe that Algorithm \ref{alg: DP-2S-GD-II} 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{sec: additional Angrist exp}.
\section{Conclusion}
We have introduced \texttt{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{rmk:compto2sls}). Improving this rate (via algorithmic modifications) and establishing lower-bounds for privacy-accuracy tradeoffs for the IVaR problem are interesting future directions.
\section*{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.
\bibliographystyle{unsrtnat}
\bibliography{references}