Optimal Algorithms for
Differentially Private  Stochastic
Monotone Variational Inequalities and Saddle-Point Problems

Digvijay Boob Engineering Management, Information, and Systems, Southern Methodist University,
Cristóbal Guzmán Department of Applied Mathematics, University of Twente. Institute for Mathematical and Computational Eng., Pontificia Universidad Católica de Chile.

In this work, we conduct the first systematic study of stochastic variational inequality (SVI) and stochastic saddle point (SSP) problems under the constraint of differential privacy (DP). We propose two algorithms: Noisy Stochastic Extragradient (NSEG) and Noisy Inexact Stochastic Proximal Point (NISPP). We show that a stochastic approximation variant of these algorithms attains risk bounds vanishing as a function of the dataset size, with respect to the strong gap function; and a sampling with replacement variant achieves optimal risk bounds with respect to a weak gap function. We also show lower bounds of the same order on weak gap function. Hence, our algorithms are optimal.

Key to our analysis is the investigation of algorithmic stability bounds, both of which are new even in the nonprivate case.

The dependence of the running time of the sampling with replacement algorithms, with respect to the dataset size $n$, is $n^2$ for NSEG and $\wt{O}(n^{3/2})$ for NISPP.

§ INTRODUCTION

Stochastic variational inequalities (SVI) and stochastic saddle-point (SSP) problems have become a central part of
the modern machine learning toolbox. The main motivation behind this line of research is
the design of algorithms for multiagent systems and adversarial training, which are more suitably modeled by
the language of games, rather than pure (stochastic) optimization. Applications that rely on these methods
may often involve the use of sensitive user data, so it becomes important to develop algorithms for these
problems with provable privacy-preserving guarantees. In this context, differential privacy (DP) has become
the gold standard of privacy-preserving
algorithms, thus a natural question is whether it is possible to design DP algorithms for SVI and SSP that attain high

Motivated by these considerations, this work provides the first systematic study of differentially-private SVI and SSP
problems. Before proceeding
to the specific results, we present more precisely the problems of interest. The stochastic variational inequality (SVI)
problem is: given a monotone
operator $F:\W\mapsto\RR^d$ in expectation form $F(w)=\EE_{\bbeta\sim {\cal P}}[F_{\bbeta}(w)]$, find
$w^* \in \W$ such that

\begin{equation}\label{main_problem}
\inprod{F(w^*)}{w-w^*} \ge 0 \quad \forall w \in \W. \tag{VI(F)}
\end{equation}

The closely related stochastic saddle point (SSP) problem is: given a convex-concave real-valued function $f:{\cal W}\mapsto\RR$
(here ${\cal W}={\cal X}\times{\cal Y}$ is a product space),
given in expectation form  $f(x,y) = \EE_{\bbeta\sim\Pcal}[f_{\bbeta}(x, y)]$, the goal is to find $(x^{\ast},y^{\ast})$
that solves

\begin{equation}\label{spp_problem}
\min_{x \in \X} \max_{y \in \Y}  f(x,y) . \tag{SP(f)}
\end{equation}

In both of these problems, the input to the algorithm is an i.i.d. sample $\bS=(\bbeta_1,\ldots,\bbeta_n)\sim {\cal P}^n$.
Uncertainty introduced by a finite random sample renders the computation of exact solutions infeasible, so gap (a.k.a. population risk)
functions are used to quantify the quality of solutions.
Let ${\cal A}:{\cal Z}^n\mapsto{\cal W}$ be an algorithm for SVI problems  (<ref>).

We define the strong VI-gap associated with ${\cal A}$ as

\begin{equation} \label{eqn:Minty_VI}
\VIgap{{\cal A}}{F}:= \added{\EE_{{\cal A},\bS}}\bracket[\Big]{\sup_{w\in\W}\langle F(w),{\cal A}(\bS)-w \rangle}.
\end{equation}

We also define the weak VI-gap as

\begin{equation} \label{eqn:weak_VI_gap}
\wVIgap{{\cal A}}{F} := \EE_{\cal A} \sup_{w \in \W} \EE_{\bS}\bracket*{\inprod{F(w)}{{\cal A}(\bS) -w}}.
\end{equation}

Here, expectation is taken over both the sample data $\bS$ and the internal randomization of ${\cal A}$.
For SSP (<ref>), given an algorithm ${\cal A}:{\cal Z}^n\mapsto {\cal X}\times{\cal Y}$, and letting ${\cal A}(\bS)=(x(\bS),y(\bS))$, a natural gap function is the following saddle-point (a.k.a. primal-dual) gap

\begin{equation} \label{eqn:SP_gap}
\SPgap{{\cal A}}{f}:= \EE_{{\cal A},\bS}\bracket[\Big]{\sup_{x\in\X,\, y\in \Y} [f(x(\bS), y)-f(x,y(\bS))] }.
%:= \EE_{{\cal A},\bS}\bracket[\Big]{\sup_{x\in\X,\, y\in \Y} [f(u_x, y)-f(x,u_y)] }.  %\tag{$\mbox{SP}(F)$}
\end{equation}

Analogously as above, we define the weak SSP gap as[The denominations of weak and strong gap functions used in this paper are not standard, but we believe are the most appropriate in this context. For example, in [41] used the terms weak and strong generalization measure for (<ref>) and (<ref>) respectively, but it is clear that these quantities do not refer to standard generalization measures used in stochastic optimization.]

\begin{equation} \label{eqn:weak_SP_gap}
\wSPgap{{\cal A}}{f}:= \EE_{\cal A} \sup_{x\in\X,\, y\in \Y} \EE_{\bS}[f(x(\bS), y)-f(x,y(\bS))] .  %\tag{$\mbox{SP}(F)$}
\end{equation}

It is easy to see that in both cases the gap is always nonnegative, and any exact solution must have zero-gap. For examples and applications of SVI and SSP we refer to Section <ref>.
Despite the fact that the strong VI is a more classical and well-studied quantity, the weak VI gap has been observed to be useful in various contexts. We refer the reader to [41] for more discussions on the weak VI gap.

On the other hand, we are interested in designing algorithms that are differentially private. These algorithms build a solution
based on a given dataset $\bS$ of random i.i.d. examples from the target distribution, and output a (randomized) feasible solution, ${\cal A}(\bS)$.
We say that two datasets $\bS=(\bbeta_i)_i,\bS^{\prime}=(\bbeta_i^{\prime})_i$ are neighbors, denoted $\bS\simeq \bS^{\prime}$, if they only differ in a single entry $i$. We say that an algorithm ${\cal A}(\bS)$ is $(\priva,\privb)$-differentially private if for every event $E$ in the output
space[Note that the probabilities in the definition of DP only involve the probability space of algorithmic randomization, and not of the datasets, which is emphasized by the notation $\PP_{\cal A}$. The datasets must be neighbors, but they are
otherwise arbitrary, and this is crucial to certify the privacy for any user.]

\begin{equation} \label{eqn:DP}
\PP_{\cal A}[{\cal A}(\bS)\in E] \leq e^{\priva}\PP_{\cal A}[{\cal A}(\bS^{\prime})\in E] + \privb
\quad(\forall \bS\simeq \bS^{\prime}).
\end{equation}

Here $\priva,\privb\geq 0$ are prescribed parameters that quantify the privacy guarantee. Designing DP algorithms for
particular data analysis problems is an active area of research. Optimal risk algorithms
for stochastic convex optimization have only very recently been developed, and it is unclear whether these methods
are extendable to SVI and SSP settings.

§.§ Summary of Contributions

Our work is the first to provide population risk bounds for DP-SVI and DP-SSP problems. Moreover, our algorithms attain
provably optimal rates and are computationally efficient. We summarize our contributions as follows:

* We provide two different algorithms for DP-SVI and DP-SSP: namely, the noisy stochastic extragradient method
() and a noisy inexact stochastic proximal-point method (NISPP). The   method is a natural DP variant
of the well-known stochastic extragradient method [22], where privacy is obtained by Gaussian noise
addition; on the other hand, the NISPP method is an approximate proximal point algorithm
[33, 20] in which every proximal iterate is made noisy to make it differentially private. Our more basic variants of both of these methods are based
on iterations involving disjoint sets of datapoints (a.k.a. single pass method), which are known to typically lead to highly suboptimal rates in DP (see the Related Work
Section for further discussion).

* We derive novel uniform stability bounds for the   and NISSP methods.
For , our stability upper bounds are inspired by the interpretation of the extragradient method as a (second order)
approximation of the proximal point algorithm. In particular, we provide expansion bounds for the extragradient iterates, and
solve a (stochastic) linear recursion. The stability bounds for NISPP method are based on stability of the (unique) SVI solution
in the  strongly monotone case. Finally, we investigate the risk attained by multipass versions of the   and NISPP methods, leveraging known generalization bounds for stable algorithms [26].

Here, we show that
the optimal risk for DP-SVI and DP-SSP can be attained by running these
algorithms with their sampling with replacement variant. In particular,  method requires $n^2$ stochastic operator evaluations, and NISPP method requires much smaller $\wt{O}(n^{3/2})$ operator evaluations for both DP-SVI and DP-SSP problems. In particular, these upper bounds also show the dependence of the running time of each of these algorithms w.r.t. the dataset size.

* Finally, we prove lower bounds on the weak gap function for any DP-SSP and DP-SVI algorithm, showing the risk optimality of the aforementioned multipass algorithms. The main challenge in these lower bounds is showing that existing constructions of lower bounds for DP convex optimization [5, 36, 4] lead to lower bounds on the weak gap of a related SP/VI problem.

The following table provides details of population risk and operator evaluation complexity .

2c|Type of sampling    2c|Type of sampling
single pass    multipass   single pass    multipass
Method   Criterion - Strong Gap    Criterion - Weak Gap    2c|Number of operator evaluations
NSEG    $O\Big(\frac{d^{1/4}}{\sqrt{n\varepsilon}}+ \frac{\sqrt{d}}{n\varepsilon} \Big)$    $O\Big(\frac{1}{\sqrt{n}} + \frac{\sqrt{d}}{n\varepsilon} \Big)$    $n$     $n^2$
NISPP (OE subroutine)    $O\Big(\frac{1}{n^{1/3}}+ \frac{\sqrt{d}}{ n^{2/3}\varepsilon}\Big)$    $O\Big(\frac{1}{\sqrt{n}} + \frac{\sqrt{d}}{n\varepsilon}\Big)$    $O(n\log{n})$    $O(n^{3/2}\log{n})$

Different levels of risk and complexity achieved  by NSEG and NISPP methods for $(\varepsilon,\eta)$-differentially private SVI/SSP. Here $n$ is the dataset size, and $d$ is the dimension of the solution search space. We omit the dependence on other problem parameters (e.g., Lipschitz constants and diameter), as well as the privacy parameter $\eta$.

§.§ Related Work

We divide our discussion on related work in three main areas. Each of these areas has been extensively investigated, so
a thorough description of existing work is not possible. We focus ourselves on the work which is more directly related to our own.

* Stochastic Variational Inequalities and Saddle-Point Problems: Variational inequalities and saddle-point problems are classical topics in applied mathematics, operations research and engineering (e.g., [32, 35, 33, 24, 30, 12, 29, 2, 31]). Their stochastic
counterparts have only gained traction recently, mainly motivated by their applications in machine learning
(e.g.,[22, 21, 17, 16, 25] and references therein).

For the stochastic version of (<ref>), [28] proposed a robust stochastic approximation method.

The first optimal algorithm for SVI with monotone Lipschitz operators was obtained by Juditsky, Nemirovski and Tauvel
[22], and very recently Kotsalis, Lan and Li [25] developed optimal variants for the strongly monotone case
(in terms of distance to the optimum criterion, rather than VI gap).

It is important to note that naive adaptation of these methods to the DP setting requires adding noise to the operator
evaluations at every iteration, which substantially degrades the accuracy of the obtained solution. A careful privacy accounting
and minibatch schedule can lead to optimal guarantees for single-pass methods [13], however
this requires accuracy guarantees for the last iterate, which is currently an open problem for SVI and SSP
(aside from specific cases, typically involving strong monotonicity conditions, e.g., [16, 25]).

We circumvent this problem by providing population risk guarantees for multipass methods.

* Stability and Generalization: Deriving generalization (or population risk) bounds for general-purpose algorithms is
a challenging task, actively studied in theoretical machine learning. Bousquet and Elisseeff [6] provided a
systematic treatment of this question for algorithms which are stable, with respect to changes of a single element
in the training dataset, and a sequence of works have refined these generalization guarantees (see [14, 7] and references therein). This idea has been applied to investigate the generalization properties of regularized empirical risk
minimization [6, 34], and more recently to iterative methods, such as stochastic gradient descent
[15, 3].

Using stability to obtain population risk bounds in SVI and SSP is substantially more challenging, due to the presence of a
supremum in the accuracy measure (see eqns. (<ref>) and (<ref>)). Recently, Zhang et al.
[41], established stability implies generalization results for the strong SP gap
under strong monotonicity assumptions.
Their proof strategy applies analogously to address the SVI setting,
although this is not carried out in their work. More recently, Lei et al. [26], proved generalization bounds on the weak
SP gap without strong monotonicity assumptions. We leverage this result for our
algorithms, and further elaborate on its implications for SVI in Section <ref>.

* Differential Privacy in (Stochastic) Convex Optimization: Differentially private empirical risk minimization and stochastic convex
optimization have been extensively studied for over a decade (see, e.g. [8, 9, 18, 23, 5, 37, 19, 4, 13]). Relevant to our work are the first optimal risk algorithms for
DP-ERM [5] and DP-SCO [4]. To the best of our knowledge, our work is the first to address
DP algorithms for SVI and SSP. Our approach for generalization of multipass algorithms is inspired by the noisy SGD analysis in [3]. However,
our stability analysis differs crucially from [3]: in the case of NSEG, we need to carefully address the double operator evaluation of the extragradient step, which is done by using the fact that the extragradient operator is approximately nonexpansive. In the case of NISPP, we leverage the contraction properties of strongly monotone VI solutions. By contrast, SGD in the nonsmooth case is far from nonexpansive [3].
Alternative approaches to obtain  optimal risk in DP-SCO, including privacy amplification by iteration [13], and
phased regularization or phased SGD [13], appear to run into fundamental limitations when applied to DP-SVI
and DP-SSP. It is an interesting future research direction to obtain faster running times with optimal population risk in DP-SVI and
DP-SSP, which
may benefit from these alternative approaches.

The main body of this paper is organized as follows. In Section <ref>, we provide the necessary background information on
SVI/SSP, uniform stability, and differential privacy, which are necessary for the rest of the paper. In Section <ref> we introduce the
method, together with its basic privacy and accuracy guarantee for a single pass version.
Section <ref> provides stability bounds for NSEG method along with the consequent optimal rates for SVI and SSP. In Section <ref>, we introduce the single-pass differentially private NISPP method with bound on expected SVI-gap. Section <ref> presents stability analysis of NISPP, together with the resulting optimal rates for
SVI/SSP gap. We conclude in Section <ref> with lower bounds that prove the optimality of the obtained rates.

§ NOTATION AND PRELIMINARIES

We work on the Euclidean space $(\RR^d, \inprod{\cdot}{\cdot})$, where $\inprod{\cdot}{\cdot}$ is

the standard inner product, and $\|u\|=\sqrt{\langle u,u\rangle}$ is the $\ell_2$-norm. Throughout, we consider a compact convex set
$\W\subseteq\RR^d$ with diameter $D>0$.  We denote the standard Euclidean projection operator on set $\W$ by $\proj_{\W}(\cdot)$.
The identity matrix on $\RR^d$ is denoted by $\mathbb{I}_d$.

We let $\mathcal{P}$ denote an unknown distribution supported on an arbitrary set $\Z$, from which we have access to exactly $n$ i.i.d. datapoints which we denote by sample set $\bS\sim {\cal P}^n$. Throughout, we will use boldface characters to denote sources of randomness (coming from the data, or internal algorithmic randomization).
We say that two datasets $\bS, \bS'$ are adjacent (or neighbors), denoted by $\bS \simeq \bS'$, if they differ in a single data point. We also denote subsets (a.k.a. batches), or single
data points,

of $\bS$ or $\Pcal$ by $\bB$ and $\bbeta$, respectively. Whether $\bbeta$ or $\bB$ is sampled from $\Pcal$ or $\bS$ is specified explicitly unless it is clear from the context. For a batch $\bB$, we denote its size by $\abs{\bB}$. Therefore, we have $\abs{\bS} = n$. Throughout, we will denote Gaussian random variables by $\boldsymbol{\xi}$.

We say that $F:\W \to \RR^d$ is a monotone operator if

\[\inprod{F(w_1) -F(w_2)}{w_1-w_2} \ge 0, \quad \forall w_1, w_2 \in \W.\]

Given $L>0$, we say that $F$ is  $L$-Lipschitz continuous, if

\[\gnorm{F(w_1) - F(w_2)}{}{} \le L\gnorm{w_1-w_2}{}{}, \quad \forall w_1, w_2 \in \W.\]

Finally, we say that $F$ is $M$-bounded if $\sup_{w\in \W}\|F(w)\|\leq M$. We denote the
set of monotone, $L$-Lipschitz and $M$-bounded operators by ${\cal M}_{\cal W}^1(M,L)$.
In this work, we will focus on the case where $F$ is an expectation operator, i.e., $F(w) := \EE_{\bbeta\sim \Pcal} [F_{\bbeta}(w)]$, where $\Pcal$ is an arbitrary distribution supported on $\Z$,

and for any $\bbeta$ in ${\cal Z}$, $F_{\bbeta}(\cdot)\in {\cal M}_{\cal W}^1(M,L)$, $\bbeta$-a.s.[
Here, we mean that for almost every $\bbeta$, we have $F_{\bbeta}\in {\cal M}_{\cal W}^1(M,L)$.

In the stochastic saddle point problem (<ref>), we modify the notation slightly. Here, ${\cal X}\subseteq\RR^{d_1}$ and ${\cal Y}\subseteq\RR^{d_2}$ are compact convex sets, and we will assume that the saddle point
functions $f_{\bbeta}(\cdot,\cdot):{\cal X}\times{\cal Y}\mapsto\RR$, satisfy the following conditions $\bbeta$-a.s.

* $\nabla_x f_{\bbeta}(\cdot,\cdot)$ is $L_x$-Lipschitz continuous, and $\nabla_y f_{\bbeta}(\cdot,\cdot)$ is $L_y$-Lipschitz continuous, and;

* $f_{\bbeta}(\cdot,y)$ is convex, for any given $y\in \Y$, and $f_{\bbeta}(x,\cdot)$ is concave, for any given $x\in \X$ (we will say in this case the function is convex-concave).

If the assumptions above are met, we will denote $L\triangleq \sqrt{L_x^2 + L_y^2}$.

It is easy to see that for any $u \in \W$, we have

$$\VIgap{u}{F}  = \sup_{w \in \W} \inprod{F(w)}{u-w} \ge \inprod{F(u)}{u-u}= 0.$$
Also, along with the above relation, monotonicity of $F$ implies that
$$ 0 \le \VIgap{w^*}{F} = \sup_{w\in \W}\inprod{F(w)}{w^*-w} \le \sup_{w\in \W} \inprod{F(w^*)}{w^*-w} \le 0.$$
The above relation yields that $\VIgap{w^*}{F} = 0$ is necessary for $w^*$ to be a solution of (<ref>). Hence, our goal in this work is to compute
a vector with the smallest possible VI gap as defined above, where the input is $\bS\sim {\cal P}^n$, and with the constraint of the algorithm being differentially private w.r.t. $\bS$ (see Section <ref> for details on DP).

For convex-concave (<ref>) problem,
our goal is to reduce the saddle-point (a.k.a. primal-dual) gap

\begin{equation} \label{eqn:SP_gap}
\SPgap{u_x}{u_y}{f}:= \sup_{x\in\X,\, y\in \Y} [f(u_x, y)-f(x,u_y)].  %\tag{$\mbox{SP}(F)$}
\end{equation}

Note that (<ref>) is slightly stronger criterion of convergence compared to (<ref>) when $F$ is the gradient operator of saddle point function $f$ (as discussed in the introduction). Indeed, when $u = \begin{bmatrix}
u_x \\ u_y
\end{bmatrix}$ is the approximate solution for associated VI-problem then convergence criterion (<ref>) implies

\begin{align*}
\sup_{w \in \W} \inprod{F(w)}{u - w} &= \sup_{x\in\X,\, y\in \Y} \inprod*{\begin{bmatrix}
\grad_x f(x, y)\\ -\grad_yf(x,y)
\end{bmatrix}}{\begin{bmatrix}
u_x -x\\ u_y-y
\end{bmatrix}}\\
&= \sup_{x\in\X,\, y\in \Y} \inprod{\grad_xf(x,y)}{u_x-x} - \inprod{\grad_yf(x,y)}{u_y-y}\\
&\le \sup_{x\in\X,\, y\in \Y} f(u_x, y) - f(x, u_y)
\end{align*}

where last inequality follows from the convex-concave nature of $f$. Not withstanding the stronger convergence criterion (<ref>), our results for (<ref>) problem can be extended to (<ref>) problem.
Under the assumptions above, it is well-known that SSP [32, 35] (and
SVI [12], respectively) have a solution.

In the case of saddle-point problems, given the convex-concave function $f_{\bbeta}(\cdot,\cdot):{\cal X}\times{\cal Y}\mapsto\RR $, it is well-known that the operator $F:{\cal X}\times{\cal Y}\mapsto \RR^d\times\RR^d$ below is monotone

\begin{equation} \label{eqn:monotone_op_SP}
F_{\bbeta}(x,y) = (\nabla_x f_{\bbeta}(x,y), -\nabla_y f(x,y)).
\end{equation}

We will call this operator the monotone operator associated with $f_{\bbeta}(\cdot,\cdot)$.
Furthermore, if $\nabla_x f_{\bbeta}(\cdot,y)$ has $L_x$-Lipschitz continuous gradient and $\nabla_y f_{\bbeta}(x,\cdot)$ has $L_y$-Lipschitz continuous gradient, then $F$ is $\sqrt{L_x^2+L_y^2}$-Lipschitz continuous.

It is easy to see that, given a SSP problem with function $f_{\bbeta}(\cdot,\cdot)$ and sets ${\cal X}$, ${\cal Y}$, an (exact) SVI solution (<ref>) for the monotone operator associated to
$f(x,y)=\EE_{\bbeta}[f_{\bbeta}(x,y)]$ over the set $\W={\cal X}\times{\cal Y}$, yields an exact SSP solution for the starting problem. Unfortunately, such reduction does not directly work for approximate solutions to (<ref>) and (<ref>), so the analysis must be done separately for both problems.

For batch $\bB$, we denote the empirical (a.k.a. sample average) operator
$F_{\bB}(w) := \frac{1}{\abs{\bB}} \sum_{\bbeta \in \bB} F_{\bbeta}(w).$
On the other hand, for a batch $\bB$, the empirical saddle point function is denoted as
$f_{\bB}(x, y) = \frac{1}{\abs{\bB}} \sum_{\bbeta \in \bB} f_{\bbeta}(x, y).$
Given a distribution ${\cal P}$, the expectation operator and function are denoted by $F_{\cal P}(w):=\mathbb{E}_{\bbeta\sim\cal P}[F_{\bbeta}(w)]$, and $f_{\cal P}(x,y)=\EE_{\bbeta\sim {\cal P}}[f_{\bbeta}(x,y)]$, respectively.
For brevity, whenever it is clear from context we will drop the dependence on ${\cal P}$.

Variational inequalities with monotone operators provide a framework for unified treatment of convex optimization and convex-concave saddle point problems.

where $\X, \Y$ are convex compact sets, $f(\cdot, y)$ is a convex function for all $y \in \Y$ and $f(x, \cdot)$ is a concave function for all $x \in \X$. Then, denoting $\W = \X \times \Y$, $w = \begin{bmatrix}
x\\ y
\end{bmatrix}$ and $F(w) = \begin{bmatrix}
\grad_x f(x,y)\\ -\grad_y f(x,y)
\end{bmatrix}$, we note that $F$ is a monotone operator and definition of $w^*$ is precisely the definition of saddle point for (<ref>) problem.

It is worth mentioning that the convex-concave function $f$ arises from an expectation $f(x,y) := \EE_{\bbeta \sim \Pcal} f(x,y;\bbeta)$, where $f(\cdot,\cdot;\bbeta)$ is convex-concave, $M$-Lipschitz and $L$-smooth, for a.e. $\bbeta$. Hence, we are essentially studying a stochastic saddle point (SSP) problem.

§.§ Examples and Applications of SVI and SSP

An interesting problem which can be formulated as a SSP-problem is the minimization of a max-type convex function:

\[\min_{x \in \X} \Big\{\phi(x) :=\max_{1 \le j \le m} \phi_j(x)\Big\},\]

where $\phi_j: \X \to \RR$ is a stochastic convex function $\phi_j(x) := \EE_{\boldsymbol{\zeta}_j\sim\Pcal_j}[\phi_{j,\boldsymbol{\zeta}_j}(x)]$ for all $j \in [m]$. This problem is essentially a structured nonsmooth optimization problem which can be reformulated into a convex-concave saddle point problem:

\[\min_{x \in \X} \max_{y \in \Delta_m} \EE_{\boldsymbol{\zeta}_1\dots\boldsymbol{\zeta}_m} [\tsum_{j=1}^m y_i\phi_{j,\boldsymbol{\zeta}_j}(x)]\]

Here, $\bbeta = (\boldsymbol{\zeta}_j)_{j=1}^m$ is the random input to the saddle point problem: $f_{\bbeta}(x,y) = \sum_{j = 1}^my_j\phi_{j,\boldsymbol{\zeta}_j}(x)$. Note that a substantial generalization of the max-type problem above is the so called compositional optimization problem:

\[\min_{x \in \X}\phi(x) := \Phi(\phi_1(x), \dots, \phi_m(x)),\]

where $\phi_j(x)$ are convex maps and $\Phi(u_1, \dots, u_m)$ is a real-valued convex function whose Fenchel-type representation
is assumed to have the form

\[\Phi(u_1, \dots, u_m) = \max_{y \in \Y} \tsum_{j=1}^m \inprod{u_j}{A_jy+ b_j} - \Phi_{*}(y),\]

where $\Phi_{\ast}$ is a convex, Lipschitz and smooth.
Then, overall optimization problem can be reformulated as a convex-concave saddle point problem:

\[\min_{x \in \X} \max_{y \in \Y} \tsum_{j=1}^m \inprod{\phi_j(x)}{A_jy+ b_j} - \Phi_{\ast}(y),\]

where stochasticity is introduced due to constituent functions $\phi_j(x) = \EE_{\boldsymbol{\zeta}_j} [\phi_{j,\boldsymbol{\zeta}_j}(x)].$

To conclude, we remark that these type of models have been recently proposed in machine learning to address approximate
fairness [39] and federated learning on heterogeneous populations [27]. In these examples, the different
indices $j\in [m]$ may denote different subgroups from a population, and we are interested in bounding the (excess) population risk
on these subgroups uniformly (with the motivation of preventing discrimination against any subgroup). This clearly cannot be achieved
by a stochastic convex program, and a stochastic saddle-point
formulation is effective in certifying accuracy across the different subgroups separately.

For further examples and applications of stochastic variational inequalities and saddle-point problems, we refer the reader to [22, 21, 41].

§.§ Algorithmic Stability

In general, an algorithm is a randomized function mapping datasets to candidate solutions, ${\cal A}:\Z^n\mapsto\RR^d$, which is measurable w.r.t. the dataset.
Two datasets, $\bS=(\bbeta_1,\ldots,\bbeta_n),\,\bS^{\prime}=(\bbeta_1^{\prime},\ldots,\bbeta_n^{\prime})\in\Z^n$ are said to be neighbors (denoted $\bS\simeq\bS^{\prime}$) if they only differ in at most one data point, namely

$$ \bbeta_j=\bbeta_j^{\prime} \qquad (\exists i\in [n])(\forall j\neq i). $$

Algorithmic stability is a notion of sensitivity analysis of an algorithm under neighboring datasets. Of particular interest to our work
is the notion of uniform argument stability (UAS).

Let ${\cal A}:\Z^n\mapsto\RR^d$ be a randomized mapping and $\delta>0$. We say that ${\cal A}$ is $\delta$-uniformly argument stable (for short, $\delta$-UAS) if

$$ \sup_{\bS\simeq\bS^{\prime}} \EE_{\cal A}\|{\cal A}(\bS)-{\cal A}(\bS^{\prime})\|\leq \delta. $$

Ocassionally, we may denote $\delta_{\cal A}(\bS,\bS^{\prime})\triangleq \|{\cal A}(\bS)-{\cal A}(\bS^{\prime})\|$, for convenience.
The importance of algorithmic stability in machine learning comes from the fact that stability implies generalization in stochastic optimization and stochastic saddle point (SSP) problems [6, 7, 41]. Below, we restate existing results on stability implies generalization for SSP problems below. Before doing so we need to briefly introduce the (strong) empirical gap function: given a dataset $\bS$ and an algorithm $\A$, we define the empirical
gap function for a saddle point and variational inequality problem respectively as

\begin{eqnarray}
\EmpSPgap{\A}{f_{\bS}} &:=& \EE_{\cal A}[\sup_{x,y} f_{\bS}(x(\bS),y)-f_{\bS}(x,y(\bS)) ]\label{eqn:emp_SP_gap}\\
\EmpVIgap{\A}{F_{\bS}} &:=& \EE_{\cal A}[\sup_{w} \langle F_{\bS}(w),\A(\bS)-w\rangle]. \label{eqn:emp_VI_gap}
\end{eqnarray}

Notice that in these definitions the dataset $\bS$ is fixed.

Consider the stochastic saddle point problem (<ref>) with functions $f_{\bbeta}(\cdot, y)$ and $f_{\bbeta}(x, \cdot)$ being $M$-Lipschitz for all $x \in\X, y \in \Y$ and $\bbeta$-a.s.. Let $\A: \Z^n \to \X \times \Y$ be an algorithm, where $\A(\bS)=(x(\bS),y(\bS))$.
If $x(\cdot)$ is $\delta_x$-UAS and $y(\cdot)$ is $\delta_y$-UAS, and both are integrable, then

\begin{equation}\label{eq:stab_gen_ssp_final}
\wSPgap{\A}{f} \le \EE_{\bS}[ \EmpSPgap{\A}{f_{\bS}} ] + M[\delta_x+\delta_y].
\end{equation}

This result can be extended for SVI problems as well. We provide a formal statement below and prove it in Appendix

Consider a stochastic variational inequality with $M$-bounded operators $F_{\bbeta}(\cdot):\W\mapsto\RR^d$.
If ${\cal A}:\Z^n\mapsto\W$ is integrable and $\delta$-UAS, then

\begin{equation}
%\EE_{{\cal A},\bS}[
\wVIgap{{\cal A}}{F}
\leq
%\EE_{{\cal A},\bS}\Big[
\EE_{\bS}[\EmpVIgap{{\cal A}}{F_\bS} ]
\end{equation}

§.§ Technical Overview

First, we propose a Noisy Stochastic ExtraGradient (NSEG) method (c.f. Algorithm <ref> and (<ref>)) for DP-SVI and DP-SSP problem. This algorithm is based on stochastic extragradient method of <cit.> where each stochastic operator update is made `noisy' by addition of Gaussian noise to achieve differential privacy. In a single pass version of this algorithm, i.e., each datapoint is used exactly once, each stochastic operator uses an independent datapoint for the evaluation. Hence, one can apply the analysis of <cit.> to obtain bound on the expected value of $\VIgap{\cdot}{F}$. For a particular step-size and batch-size policy, we show that NSEG method achieves a bound of $O(\max\{\frac{(d\log(1/\delta))^{1/4}}{\sqrt{n}}, \frac{\sqrt{d\log(1/\delta)}}{\theta n} \})$ on VI-gap and SP-gap, where $(\theta, \delta)$ are constants associated with the privacy budget. We also provide a different step-size and batch-size policy for which NSEG method achieves $O(\frac{1}{n^{1/3}} + \frac{\sqrt{d \log(1/\delta)}}{\theta n^{2/3}})$ bound on the VI-gap and SP-gap. This suggest that if we can run the algorithm for $O(n^{3/2})$ stochastic operator evaluations (much longer than $n$ stochastic operator evaluations in single pass regime) then there is a possibility that we can achieve optimal convergence guarantee.

However, when running this algorithm for longer iterations, we loose the fact that each stochastic operator has an independently sampled datapoint from distribution $\Pcal$. We call this a  multipass regime over dataset $\bS$. A direct application of previous results can only guarantee a bound on expected VI-gap of sample average operator, i.e., a bound $\VIgap{\cdot}{F(\cdot; \bS)}$. In order to convert this guarantee into $\VIgap{\cdot}{F}$, we perform generalization analysis of this algorithm using stability arguments.

“Stability implies generalization" results are well-known for convex optimization. In particular, [6] showed such results for deterministic optimization of the sample average problem and [15] for stochastic optimization of the sample average problem. Such results are difficult to extend to SSP or SVI problems due to presence of $\sup$ in both convergence criteria (<ref>) and (<ref>). To our best knowledge, the only result in the literature for generalization of SSP is in [41] where there results are for exact solution of the subproblem and the criterion is significantly weaker than the one in  (<ref>), i.e., our criterion is

\begin{align*}
\EE_{{\cal A},\bS} [\SPgap{u_x(\A, \bS)}{u_y(\A, \bS)}{f}] = \EE_{{\cal A},\bS} [\sup_{x\in\X,\, y\in \Y} f(u_x(\A, \bS), y) - f(x, u_y(\A, \bS))]
\end{align*}

whereas in [41] assumes that algorithm $\A$ is exact for the sample average problem (hence, can be ignored)  and take the expectation w.r.t. $\bS$ inside the supremum, i.e.,

\begin{align*}
\text{Weak}\SPgap{u_x(\bS)}{u_y(\bS)}{f} = \sup_{x\in\X,\, y\in \Y} \EE_{\bS} [f(u_x(\bS), y) - f(x, u_y(\bS))].
\end{align*}

discussion of Ozdaglar below
There is also recent work of <cit.> which discusses generalization guarantee for SSP problem. However, their current draft seems to have an error associated with handling max term (see proof of Theorem 1 in their paper).
Our analysis of “stability implies generalization” provides a novel and simple way to handle the max term and hence provides the strongest known generalization guarantee for SSP problem and first guarantee for SVI problem.

Using stability implies generalization bounds for NSEG, we provide bounds on the VI-gap and SP-gap of $O(\frac{1}{\sqrt{n}}+ \frac{\sqrt{d\log(1/\delta)}}{n})$ which is shown to be optimal risk. However, those bounds are achieved in $O(n^2)$ stochastic operator evaluations as compared to earlier estimate of $O(n^3/2)$ stochastic operator evaluations. The key reason for this additional cost is that extragradient method is not very stable. may be some citation needed in next statement. Otherwise can remove. Actually,
the bound does not get stuck because of stability, but rather from the convergence guarantee of SEG. I mean, stability also seems to be on the verge of not working, but this is because of noise addition, not because of the stability of SEG. Anyway, I think this claim is not accurate. In fact, it is well-known that accelerated first-order methods are less stable compared to their non-accelerated peers (<cit.>)

In order to get around this difficulty, we propose, to our best knowledge, the first  inexact proximal point method for VI problems. We call this Noisy Inexact Stochastic Proximal Point (NISPP) method due to addition of Gaussian noise in between phases to ensure differential privacy. For single pass version of NISPP, we show results similar to that of NSEG with additional $\log{n}$ in running time complexity. Also, we show that this method is much more stable compared to NSEG and hence can achieve the optimal bounds on VI-gap and SP-gap in $O(n^{3/2}\log{n})$ stochastic operator evaluations.

§.§ Background on differential privacy

Differential privacy is an algorithmic stability type of guarantee for randomized algorithms, that certifies that the output distribution
of the algorithm “does not change too much” by changes in a single element from the dataset. The formal definition is
provided in eqn. (<ref>).

Next we provide some basic results in differential privacy, which we will need for our work. For further information on the topic,
we refer the reader to the monograph [10].
We will also consider an alternative notion of privacy known as Rényi differential privacy <cit.>. For this, we recall the definition of the Rényi divergence of order $\alpha>1$ between two probability measures $\PP$ and $\QQ$, $D_{\RDPa}(\PP||\QQ) =\frac{1}{\RDPa-1}\log \EE_{\xi\sim\QQ}\big(\frac{d\PP}{d\QQ}\big)^{\RDPa}$, where $d\PP/d\QQ$ is the Radon-Nykodim derivative. We say that an algorithm ${\cal A}:{\cal Z}^n\mapsto {\cal X}$ is $(\RDPa,\RDPb)$-Rényi differentially private (RDP) if

$$ D_{\RDPa}({\cal A}(S)||{\cal A}(S^{\prime})) \leq \RDPb \quad (\forall S\simeq S^{\prime}).$$

An important property of the Rényi divergence is its joint quasi-convexity: if $\lambda\geq 0$, $\sum_{i\in[n]}\lambda_i =1$, and $\PP_1,\QQ_1,\ldots,\PP_n,\QQ_n$ are probability measures, then $D_{\RDPa}(\sum_{i\in[n]}\lambda_i\PP_i || \sum_{i\in[n]}\lambda_i\QQ_i) \leq \max_{i\in[n]} D_{\RDPa}(\PP_i ||\QQ_i)$ <cit.>.

Finally, RDP and DP can be related in the following way.

If ${\cal A}$ is $(\RDPa,\RDPb)$-RDP, then it is also $(\RDPb+\frac{\log 1/\delta}{\RDPa-1},\delta)$-DP, for any $0<\delta<1$.

§.§.§ Basic privacy guarantees

In this work, most of our privacy guarantees will be obtained by the well-known Gaussian mechanism,
which performs Gaussian noise addition on a function with bounded sensitivity. Given a function
${\cal A}:\Z^n\mapsto \RR^d$, we define its $\ell_2$-sensitivity as

\begin{equation} \label{eqn:sensitivity}
\sup_{\bS\simeq \bS^{\prime}}\|{\cal A}(\bS)-{\cal A}(\bS^{\prime})\|.
\end{equation}

If ${\cal A}$ is randomized, then the supremum must hold with high-probability over the
randomization of ${\cal A}$ (this will not be a problem in this work, since our randomized algorithms enjoy sensitivity bounds w.p. 1).
The Gaussian mechanism (associated to function ${\cal A}$) is defined as
${\cal A}_{\cal G}(\bS)\sim {\cal N}({\cal A}(\bS),\sigma^2I)$.

Let ${\cal A}:\Z^n\mapsto \RR^d$ be a function with $\ell_2$-sensitivity $s>0$. Then, for
$\sigma^2=2s^2\ln(1/\privb)/\priva^2$,  the Gaussian mechanism is $(\priva,\privb)$-DP.

Our algorithms will adaptively use a DP mechanism such as the above. Certifying privacy of a composition can be
achieved in different ways. The most basic result establishes that if we use disjoint batches of data at each iteration, then the
composition will preserve the largest privacy parameter among its building blocks. This result is known as
parallel composition, and its proof its a direct application of
the post-processing property of DP.

Let $\bS=(\bS_1,\ldots,\bS_K)\in\Z^n$ be a dataset partitioned on blocks of sizes $n_1,\ldots,n_K$, respectively.
${\cal A}_k:\Z^{n_k}\times \RR^{d\times (k-1)}\mapsto \RR^d$, $k=1,\ldots,K$, be a sequence of mechanisms,
and let ${\cal A}:\Z^n\mapsto\RR^d$
be given by

\begin{eqnarray*}
{\cal B}_1(\bS) &=& {\cal A}_1(\bS_1)\\
{\cal B}_{k}(\bS) &=& {\cal A}_{k}(\bS_{k}, {\cal B}_{1}(\bS),{\cal B}_{2}(\bS),\ldots,{\cal B}_{k-1}(\bS)) \quad (\forall k=2,\ldots,K-1)\\
{\cal A}(\bS)            &=&{\cal A}_{K}(\bS_{K},{\cal B}_{1}(\bS),{\cal B}_{2}(\bS),\ldots,{\cal B}_{K-1}(\bS)).
\end{eqnarray*}

If each ${\cal A}_k$ is $(\priva_k,\privb_k)$-DP in its first argument (i.e., w.r.t. $\bS_k$) then ${\cal A}$ is $(\max_k \priva_k, \max_k \privb_k)$-DP.

Some of the algorithms we develop in this work make repeated use of the data, and certifying privacy for these
algorithms requires the use of adaptive composition results in DP (see, e.g. [11, 10]). For our algorithms, it
is particularly important to leverage the sampling with replacement procedure to select the data that is used
at each iteration, for which sharp bounds on DP can be obtained by the moments accountant method [1].
Below we summarize a specific version of this method that suffices for our purposes.[In our case we use uniform sampling on each iteration, as opposed to the Poisson sampling of [1];
however, it is possible to verify that similar moment estimates lead to our stated result.]

We consider an sequence of functions ${\cal A}_1,\ldots,{\cal A}_K$, where ${\cal A}_k:\Z^{n_k}\times \RR^{d\times(k-1)}\mapsto\RR^d$
is a private mechanism w.r.t. its first input (i.e., $S_k\in \Z^{n_k}$). Next we compose these mechanisms by subsampling from the
dataset, i.e., letting  $S_k\sim(\mbox{Unif}[n])^{n_k}$ for $k=1,\ldots,K$ independently, and

\begin{eqnarray*}
{\cal B}_1(S) &=& {\cal A}_1(S_1) \\
{\cal B}_{k}(S) &=& {\cal A}_{k}(S_{k}, {\cal B}_{1}(S),{\cal B}_{2}(S),\ldots,{\cal B}_{k-1}(S)) \quad (\forall k=2,\ldots,K).
%{\cal A}(S)            &=&{\cal A}_{K}(S_{K},{\cal B}_{1}(S),{\cal B}_{2}(S),\ldots,{\cal B}_{K-1}(S)).
\end{eqnarray*}

Let $w_k={\cal B}_k(S)$, according to the recursion above, and let

$\A(\bS) := \{w_1, \dots, w_K\}$ be the complete transcript of the procedure,

where we also denote output of $i$-th iterate of Algorithm <ref> as $\A_i(\bS)$.
Then privacy loss random variable can be defined as

\begin{align*}
\boldsymbol{\Upsilon}(\bS, \bS') &:= \log{\tfrac{d\Prob[\A(\bS) =  \{w_1, \dots, w_K\}]}{d\Prob[\A(\bS') =  \{w_1, \dots, w_K\}]}} = \log {\prod_{i = 1}^K \frac{d\Prob[\A_i(\bS) = w_i| \{w_{i-1}, \dots, w_1, w_0\}]}{d\Prob[\A_i(\bS') = w_i| \{w_{i-1}, \dots, w_1, w_0\}]}},%\\
%	\boldsymbol{\Upsilon} &= \max_{\bS \sim \bS'} \boldsymbol{\Upsilon}(\bS, \bS')
\end{align*}

where the second equality follows by multiplication rule. In order to prove $(\theta, \delta)$-differential privacy, it suffices to show that for all neighboring datasets $\bS\sim\bS'$,  we have $\boldsymbol{\Upsilon} \le \theta$ w.p. at least $1 -\delta$.  In the  moment accountant technique, we consider the following moment generating function of $\boldsymbol{\Upsilon}$:

\[ \Prob[\boldsymbol{\Upsilon}(\bS, \bS') > \theta] = \Prob[e^{q \boldsymbol{\Upsilon} } > e^{q\theta}] \le \frac{\EE[\exp\{q\boldsymbol{\Upsilon}(\bS, \bS')\}]}{\exp\{q\theta\}}, \]

for any $q> 0$ and $\bS \sim \bS'$. It suffices to show that

\begin{equation}\label{eq:int_rel_SPP_7}
\inf_{q > 0}\sup_{\bS\sim \bS'}\frac{\EE[\exp\{q\boldsymbol{\Upsilon}(\bS, \bS')\}]}{\exp\{q\theta\}} \le \delta.
\end{equation}

\begin{align*}
c_i(q; \bS, \bS', \{w_{i-1}, \dots, w_0\}) &:= \log\EE_{w_i} \bracket*{\paran*{\frac{d\Prob[\A_i(\bS) = w_i| \{w_{i-1}, \dots, w_1, w_0\}]}{d\Prob[\A_i(\bS') = w_i| \{w_{i-1}, \dots, w_1, w_0\}]}}^q},\\
c_i(q; \bS, \bS') &:= \sup_{\{w_{i-1}, \dots, w_0\} \in \W^i} c_i(q; \bS, \bS', \{w_{i-1}, \dots, w_0\}),\\
c(q; \bS, \bS') &:= \tsum_{i = 1}^k c_i(q; \bS, \bS'),\\
c(q) := \sup_{\bS\sim \bS'} c(q; \bS, \bS'), &\quad\quad c_i(q) := \sup_{\bS\sim \bS'} c_i(q; \bS, \bS')
\end{align*}

By sequential conditioning, we get that

\begin{align*}
\EE[\exp\{q\boldsymbol{\Upsilon}\} ] &= \EE_{w_1, \dots, w_{K-1}}\bracket*{ \prod_{i = 1}^{K-1} \bracket*{\frac{d\Prob[\A_i(\bS) = w_i| \{w_{i-1}, \dots, w_1, w_0\}]}{d\Prob[\A_i(\bS') = w_i| \{w_{i-1}, \dots, w_1, w_0\}]} }^q \EE_{w_K}\bracket*{\frac{d\Prob[\A_i(\bS) = w_K| \{w_{K-1}, \dots, w_1, w_0\}]}{d\Prob[\A_i(\bS') = w_K| \{w_{K-1}, \dots, w_1, w_0\}]}}^q}\\
&= \EE_{w_1, \dots, w_{K-1}}\bracket*{ \prod_{i = 1}^{K-1} \bracket*{\frac{d\Prob[\A_i(\bS) = w_i| \{w_{i-1}, \dots, w_1, w_0\}]}{d\Prob[\A_i(\bS') = w_i| \{w_{i-1}, \dots, w_1, w_0\}]} }^q \exp\braces*{c_K(q; \bS, \bS', \{w_{K-1}, \dots, w_0\})}}\\
&\le \exp\braces*{c_K(q; \bS, \bS')}\EE_{w_1, \dots, w_{K-1}}\bracket*{ \prod_{i = 1}^{K-1} \bracket*{\frac{d\Prob[\A_i(\bS) = w_i| \{w_{i-1}, \dots, w_1, w_0\}]}{d\Prob[\A_i(\bS') = w_i| \{w_{i-1}, \dots, w_1, w_0\}]} }^q}\\
&\le \exp\braces*{\tsum_{i = 1}^K c_i(q; \bS, \bS')}\\
&= \exp\braces*{ c(q; \bS, \bS')},
\end{align*}

from which we conclude that

$$ \sup_{\bS\sim \bS'} \EE[\exp\{ q\boldsymbol{\Upsilon}(\bS, \bS')\}] \le \exp\{c(q)\} \le \exp\{\tsum_{i = 1}^K c_i(q)\}.$$

Hence, it suffices to bound individual $c_i(q)$ to obtain a working bound on privacy loss. Note that for Algorithm <ref>, the quantity $c_i(q; \bS, \bS', \{w_{i-1}, \dots, w_0\})$ dependence on the past iterates reduces to only $w_{i-1}$. Let us define mapping $M_i(\cdot; \bS): w_{i-1} \to \wt{w}_i$ on dataset $\bS$. Then, we have $w_i = M_i(w_{i-1}; \bS) + \boldsymbol{\xi}_i$.
Then, using standard arguments, we have privacy loss variable (or ratio of distributions $M_i(w_{i-1}; \bS) + \boldsymbol{\xi}_i$ and $M_i(w_{i-1}; \bS') + \boldsymbol{\xi}_i$) conditioned on
$w_{i-1}$ is equal to

\[ \frac{d\Prob[\A_i(\bS) = w_i|w_{i-1}]}{d\Prob[\A_i(\bS') = w_i|w_{i-1}]} = \frac{d\Prob[\boldsymbol{\xi}_i = u]}{d\Prob[\boldsymbol{\xi}_i = u + \boldsymbol{\Delta}_i]},\]

where $\boldsymbol{\Delta}_i = \gnorm{M_i(w_{i-1}; \bS) - M_i(w_{i-1}; \bS')}{}{}$. Also, note that spherically symmetric normal random variable $\boldsymbol{\xi}_i$ and any arbitrary orthogonal matrix $U \in \RR^{d\times d}$, we have $U\boldsymbol{\xi}_i$ has the same distribution as $\boldsymbol{\xi}_i$. Assume, without loss generality, that first row of $U$ is aligned with $\boldsymbol{\Delta}_i$. Then, we have

\[\frac{d\Prob[\boldsymbol{\xi}_i = u]}{d\Prob[\boldsymbol{\xi}_i = u + \boldsymbol{\Delta}_i]} = \frac{d\Prob[U\boldsymbol{\xi}_i = Uu]}{d\Prob[U\boldsymbol{\xi}_i = Uu + U\boldsymbol{\Delta}_i]} = \frac{d\Prob[\boldsymbol{\xi}_i = Uu]}{d\Prob[\boldsymbol{\xi}_i = Uu + \gnorm{\boldsymbol{\Delta}_i}{}{}e_1 ] } = \frac{d\Prob[\boldsymbol{\xi}_{i,1} = u_1]}{d\Prob[\boldsymbol{\xi}_{i,1} = u_1 + \gnorm{\boldsymbol{\Delta}_i}{}{}]} %= \exp\{-\tfrac{1}{2\sigma_i^2}[\gnorm{u}{}{2} - \gnorm{}{}{2}]\}
\]

Since there is a change in only first coordinate, the problem reduces to its one-dimensional marginal.

Moreover, we know from Proposition <ref> that
$\gnorm{\boldsymbol{\Delta}_i}{}{} \le \tfrac{4M}{B_{k+1}\lambda_k}$. Let us define $p_i = \tfrac{\abs{\bB_i}}{n}$ then,

\[ \|\boldsymbol{\Delta}_i\|  \begin{cases}
=0 &\text{ if }\bB_i = \bB_i'\\
\le \tfrac{4M}{B_{i}\lambda_{i-1}} &\text{otherwise}
\end{cases}
\quad \Longrightarrow \quad
\|\boldsymbol{\Delta}_i\| \begin{cases}
= 0 &\text{w.p. } 1 - p_i\\
\le \tfrac{4M}{B_i\lambda_{i-1}} &\text{w.p. } p_i.
\end{cases}\]

By the law of total probability, we note that the denominator of our density ratio has a mixture form:

\begin{align*}
d\Prob[\boldsymbol{\xi}_{i,1} = u_1 + \gnorm{\boldsymbol{\Delta}_i}{}{}] %&= \Prob\big[\boldsymbol{\Delta}_i = 0\big] \Prob\big[\boldsymbol{\xi}_{i,1} = u_1 + \gnorm{\boldsymbol{\Delta}_i}{}{}| \boldsymbol{\Delta}_i = 0\big]   + \Prob\big[\boldsymbol{\Delta}_i \ne 0\big] \Prob\big[\boldsymbol{\xi}_{i,1} = u_1 + \gnorm{\boldsymbol{\Delta}_i}{}{}| \boldsymbol{\Delta}_i \ne 0\big]\\
&= p_id\Prob[\boldsymbol{\xi}_{i,1} = u_1 + \gnorm{\boldsymbol{\Delta}_{i}}{}{}] + (1-p_i) d\Prob[\boldsymbol{\xi}_{i,1} = u_1 ].
\end{align*}

Hence, following the technical details in the proof of Lemma 3 in [1] (which provide bounds on moment generating function of privacy loss random variable of 1-D Gaussian mechanism) and noting that  Do we need the next equation?

\begin{equation}\label{eq:int_rel_SPP_2}
\sigma_i = \tfrac{4M}{B_i\lambda_{i-1}}\sigma,
\end{equation}

we have

\begin{align*}
c_i(q; \bS, \bS', \{w_{i-1}, \dots, w_0\}) &= \log \EE \bracket*{\paran*{ \frac{d\Prob[\A_i(\bS) = w_i|w_{i-1}]}{d\Prob[\A_i(\bS') = w_i|w_{i-1}]} }^q }
= \log \EE \bracket*{ \paran*{ \frac{d\Prob[\boldsymbol{\xi}_i = u]}{d\Prob[\boldsymbol{\xi}_i = u + \boldsymbol{\Delta}_i]} }^q }\\
&= \log \EE \bracket*{ \paran*{ \frac{d\Prob[\boldsymbol{\xi}_{i,1} = u_1]}{d\Prob[\boldsymbol{\xi}_{i,1} = u_1 + \gnorm{\boldsymbol{\Delta}_i}{}{}]} }^q }\\
&= \log \EE \bracket*{ \paran*{ \frac{\Prob[\boldsymbol{\xi}_{i,1} = u_1]}{ p_i\Prob[\boldsymbol{\xi}_{i,1} = u_1 + \gnorm{\boldsymbol{\Delta}_{i}}{}{}] + (1-p_i) \Prob[\boldsymbol{\xi}_{i,1} = u_1 ] } }^q }\\
&\le %\frac{q(q+1)p_i^2}{2(1-p_i)}\paran*{\exp\braces*{\paran*{\frac{4M}{B_i\lambda_{i-1}\sigma_i } }^2 } - 1}\\
\frac{q(q+1)p_i^2}{2(1-p_i)}\paran*{\exp\braces*{1/\sigma^2}  - 1} + O(p_i^3q^3/\sigma^3)\\
&\le \frac{q^2p_i^2}{\sigma^2}
\end{align*}

where in the first inequality we used Lemma <ref>, and

where in the last inequality, we assumed that $p_iq/\sigma \ll 1$, $q \gg 1$ and $p_i \ll 1$ be such that the bound holds. (This will be satisfied for large datasets, i.e., $n \gg 1$). Note that this bound is independent of $w_{i-1}$ as well as maximum with respect to neighboring datasets $\bS, \bS'$ since we know from Proposition <ref> that $\gnorm{\boldsymbol{\Delta}_i}{}{} \le \tfrac{4M}{B_i\lambda_{i-1}}$. Hence, we conclude that

$$c_i(q) = \max_{\bS\sim \bS'}\max_{\{w_{i-1}, \dots, w_0\} \in \W^i} c_i(q; \bS, \bS', \{w_{i-1},\dots, w_0\}) \le \frac{q^2p_i^2}{\sigma^2}.$$

Finally, noting that

\[\inf_{q > 0} \max_{\bS\sim \bS'} \EE[\exp\{q\boldsymbol{\Upsilon}(\bS, \bS') - q \theta\}] \le \inf_{q > 0} \exp\{\tsum_{i = 1}^k c_i(q) - q\theta\} \le \inf\limits_{q > 0}\exp\{ \tfrac{q^2}{\sigma^2}\tsum_{i = 1}^Kp_i^2 - q \theta\}\]

\begin{equation}\label{eq:int_rel_SPP_1}
q = \tfrac{\theta \sigma^2}{2 \tsum_{i = 1}^Kp_i^2},
\end{equation}

we have
$\inf_{q > 0} \max_{\bS\sim \bS'} \EE[\exp\{q\boldsymbol{\Upsilon}(\bS, \bS') - q \theta\}] \le \exp\{\tfrac{-q\theta}{2}\} $. Now to satisfy (<ref>), we require $\exp\{-q\theta/2\} \le \delta$ implying that $q = \frac{2}{\theta}\log{\tfrac{1}{\delta}}$. Using this value in  (<ref>), we have

\begin{equation}\label{eq:int_rel_SPP_3}
\sigma = \frac{2}{\theta}\sqrt{\log{\frac{1}{\delta}}\tsum_{i = 1}^K p_i^2}.
\end{equation}

Hence, we conclude the proof.

Let $f:\Z\mapsto\RR$ be a function with $|f(z)|\leq 1$ for all $z\in\Z$. Let $\sigma\geq 1$, $0\leq p\leq 1$ with
$p\leq \sigma^2\ln(1/[p\sigma])$. Consider the distributions $\mu_1$ and $\mu_2$ induced by the
random variables

$$ {\cal N}(0,\sigma^1) \qquad \mbox{ and } \qquad q \cdot {\cal N}(f(\Delta),\sigma^1)+(1-q)\cdot{\cal N}(0,\sigma^1), $$

respectively (where $\Delta$ is any random variable supported in $\Z$). Then

$$ \max\Big\{ \EE_{\mu_2}\Big[ \Big( \frac{d\mu_1}{d\mu_2} \Big)^q \Big], \EE_{\mu_1}\Big[ \Big( \frac{d\mu_2}{d\mu_1} \Big)^q \Big] \Big\}\leq \frac{p^2q(q+1)}{(1-q)\sigma^2}+O\Big(\frac{p^3q^3}{\sigma^3} \Big). $$

Consider sequence of functions ${\cal A}_1,\ldots,{\cal A}_K$,
where ${\cal A}_k:\Z^{n_k}\times \RR^{d\times(k-1)}\mapsto\RR^d$ is a function with
sensitivity bounded as a function of the last data batch size, as follows

$$ \sup_{L\in\RR^{d\times(k-1)},\bS_k\simeq\bS_k^{\prime}}\|{\cal A}_k(\bS_k, L)-{\cal A}_k(\bS_k^{\prime}, L)\|\leq {s}%{|\bS_k|}

Consider the mechanism obtained by sampling a random subset of size $m$ from the dataset, i.e., letting $\bS_k\sim (\mbox{Unif}([\bS]))^m$, and composing it with a Gaussian mechanism with noise $\sigma^2$, i.e.

\begin{eqnarray*}
{\cal B}_1(\bS) &=& ({\cal A}_1)_{\cal G}(\bS_{1}) \\
{\cal B}_{k}(\bS) &=& ({\cal A}_{k})_{\cal G}(\bS_{k}, {\cal B}_{1}(\bS),{\cal B}_{2}(\bS),\ldots,{\cal B}_{k-1}(\bS)) \quad (\forall k=2,\ldots,K).
%{\cal A}(S)            &=&{\cal A}_{K}(S_{K},{\cal B}_{1}(S),{\cal B}_{2}(S),\ldots,{\cal B}_{K-1}(S)).
\end{eqnarray*}

There exists an absolute constant $c_1>0$, such that if $\priva<c_1K(m/n)^2$ and the noise parameter $\sigma\geq \sqrt{2K\ln(1/\privb)}sm/[n\priva]$, then  $\A(\bS) := \{{\cal B}_1(\bS), \ldots, {\cal B}_K(\bS)\}$  is $(\priva,\privb)$-differentially private.

§ THE NOISY STOCHASTIC EXTRAGRADIENT METHOD

To solve the DP-SVI problem we propose
a noisy stochastic extragradient method (NSEG) in Algorithm <ref>.
Noisy Stochastic Extragradient (NSEG) Method

[1]
Input: Starting point $u_0 \in \W$, dataset $\bS=(\bbeta_i)_{i\in[n]}\sim{\cal P}^n$,

stepsizes $(\gamma_t)_{t\in[T]}$
$t = 1, \dots, T$

$F_{1,t}(\cdot) = F_{\bB_t^1}(\cdot)+\boldsymbol{\xi}_t^1$, where $\bB_t^1\subseteq\bS$ and  $\boldsymbol{\xi}_t^1\sim {\cal N}(0,\sigma_t^2)$
$w_t = \proj_{\W}(u_{t-1}-\gamma_t F_{1,t}(u_{t-1}))$
$F_{2,t}(\cdot)= F_{\bB_t^2}(\cdot)+\boldsymbol{\xi}_t^2$, where $\bB_t^2\subseteq\bS$ and $\boldsymbol{\xi}_t^1\sim {\cal N}(0,\sigma_t^2)$
$u_t = \proj_{\W}(u_{t-1}-\gamma_t F_{2,t}(w_t))$

$\wb{w}^T =(\tsum_{t = 1}^T\gamma_t)^{-1} \tsum_{t = 1}^{T}\gamma_tw_t$

The name noisy and stochastic in Algorithm <ref> is justified by the sequence of operators $F_{1,t},F_{2,t}$ we use:

\begin{eqnarray} \label{eqn:stoch_oracle_nseg}
F_{1,t}(\cdot) \triangleq F_{\bB_t^1}(\cdot)+\boldsymbol{\xi}_t^1,
&\quad& F_{2,t}(\cdot)\triangleq F_{\bB_t^2}(\cdot)+\boldsymbol{\xi}_t^2.
\end{eqnarray}

where $\bB_t^1,\bB_t^2$ are batches extracted from dataset $\bS$, and $\boldsymbol{\xi}_t^1,\boldsymbol{\xi}_t^2\stackrel{i.i.d.}{\sim}{\cal N}(0,\sigma_t^2)$. We will denote the batch size of batch $\bB_t^1$ and $\bB_t^2$ as

The exact details of the sampling method for $\bB_t$ will depend on the variant of the algorithm.

Here, we detail some key features of the above algorithm. Stochastic extragradient was proposed in [22] where they do not have any noise addition in $F_{1,t}, F_{2,t}$ (stochasticity only arises from the dataset randomness), and where disjoint batches
are used for all iterations, as well as within iterations. This choice
is motivated by the goal of extracting population risk bounds for their algorithm.

Another important consideration is that this algorithm can also be applied to an SSP problem by using as stochastic oracle the monotone operator associated to the stochastic convex-concave function (<ref>), over the set $\W={\cal X}\times{\cal Y}$. From here onwards, when we say that a certain SVI algorithm is applied to an SSP, we mean using the choices above for the operator and feasible set, respectively.

We start by stating the convergence guarantees for the  single-pass   method. This is obtained as a direct corollary of
<cit.>, where we use an explicit bound on the oracle error with the variance of the Gaussian.

Consider a stochastic variational inequality (<ref>), with operators $F_{\bbeta}$
in ${\cal M}^1(L,M)$. Let $\A$ be the  method (Algorithm <ref>) where $0<\gamma_t\leq 1/[\sqrt3 L]$ and $(\bB_t^1,\bB_t^2)_{t}$ are
independent random variables from a product distribution $\bB_t^1,\bB_t^2\sim{\cal P}^{B_t}$, satisfies

\begin{equation*}
%\EE_{\bS,{\cal A}}\Big[
\VIgap{{\cal A}}{F}
\leq \frac{K_0(T)}{\Gamma_T},
\end{equation*}

where $K_0(T)\triangleq\Big( D^2 +7\sum_{t\in[T]}\gamma_t^2[M^2/2+d\sigma_t^2]\Big)$,
$\Gamma_T = \tsum_{t=1}^T\gamma_t$, ${\cal A}(\bS)$ is the output of Algorithm <ref> on the dataset $\bS=\bigcup_t \bB_t \sim{\cal P}^n$ and expectation in the left hand side is taken over the dataset draws, random sample batch choices, as well as noise $\bg^t_1, \bg^t_2$.

On the other hand, $\A$ applied to a stochastic (<ref>) problem attains saddle point gap

\begin{equation*}
%\EE_{\bS,{\cal A}}\Big[
\SPgap{{\cal A}}{f}
\leq \frac{K_0(T)}{\Gamma_T}.
\end{equation*}

Our first observation states that two different forms of perturbation that one might consider to ensure privacy –namely, output perturbation
(after evaluating each update) and perturbing the operator (before computing each update)–
are nearly equivalent. This near equivalence becomes exact equivalence for the resolvent update, which corresponds to $s=\infty$.

Let $F\in{\cal M}_{\W}^1(M,L)$, $g\in \bE$ (which we interpret as a linear operator). Then

$$ \|\T{s}{\gamma(F+g)}{u}-\T{s}{\gamma F}{u-\gamma g}\| \leq 2D(\gamma L)^s. $$

We first observe that $\T{}{\gamma(F+g)}{u}=\T{}{\gamma F}{u-\gamma g}$. Indeed, they both respectively satisfy the fixed point equation $v = \proj_{\W}(u-\gamma g-\gamma F(v)),$
which has a unique solution $v$ by the uniqueness of solutions to (<ref>).

Using this fact, together with eqn. (<ref>),

\begin{eqnarray*}
&&\|\T{s}{\gamma(F+g)}{u}-\T{s}{\gamma F}{u-\gamma g}\| \\
&\leq& \|\T{s}{F+\gamma g}{u}-\T{}{\gamma(F+g)}{u}\|%+\|T(u;\gamma(F+g))-T(u-\gamma g;\gamma F)\|%\\
+\|\T{}{\gamma F}{u-\gamma g}-\T{s}{\gamma F}{u-\gamma g}\|\\
&\leq& 2D(\gamma L)^s.
\end{eqnarray*}

Now we focus on convergence properties of Algorithm <ref>. In the following lemma, we give a simple extension of a key result in [29].

Suppose $\gamma_tL \le \tfrac{1}{\sqrt{2}}$. Then, for every $t=0,\ldots,T-1$ of Algorithm <ref>, we have

$$ \langle \gamma_t F_{2,t}(w_t),w_t-u \rangle
\leq \frac12\|u-u_t\|^2-\frac12\|u-u_{t+1}\|^2+\gamma_t^2\gnorm{v_t}{}{2} \quad(\forall u\in\W),$$

where $v_t :=F_{2,t}(w_t) - F_{1,t}(w_t)$.

By the optimality conditions defining $u_{t}$, we have for all $u\in\W$

\begin{eqnarray}
&\langle \gamma_tF_{2,t}(w_t),u_{t}-u\rangle
&\leq \tfrac{1}{2}\bracket*{\gnorm{u-u_{t-1}}{}{2} - \gnorm{u-u_{t}}{}{2}-\gnorm{u_{t-1}-u_t}{}{2}} \notag\\
\Rightarrow & \langle \gamma_tF_{2,t}(w_t),w_t-u\rangle
&\le\tfrac{1}{2}\bracket*{\gnorm{u-u_{t-1}}{}{2} - \gnorm{u-u_{t}}{}{2}-\gnorm{u_{t-1}-u_t}{}{2}}
\label{eqn:last_step_EG_new}
\end{eqnarray}

Next, by the optimality conditions that define $w_t$, we have

\begin{align*}
\langle \gamma_t F_{1,t}(u_{t-1}),w_t-u\rangle &\leq
\tfrac{1}{2}\bracket*{ \gnorm{u - u_{t-1}}{}{2}-\gnorm{u -w_t}{}{2}-\gnorm{w_t-u_{t-1}}{}{2}}. %\\
%\langle \gamma_t F_{1,t}(u_t),u_1^t-u_{t+1}\rangle &\leq& D(u_{t+1},u_t)-D(u_{t+1},u_1^t)-D(u_1^t,u_t).
\end{align*}

Choosing $u = u_t$ in the above relation and then summing with (<ref>), we obtain

\begin{align*}
\gamma_t \inprod{F_{2,t}(w_t)}{w_t-u} &+ \gamma_t\inprod{F_{2,t}(w_t)-F_{1,t}(u_{t-1})}{u_t-w_t} \\
&\quad\le \tfrac{1}{2}\bracket*{ \gnorm{u-u_{t-1}}{}{2} - \gnorm{u-u_{t}}{}{2} -\gnorm{u_t -w_t}{}{2}-\gnorm{w_t-u_{t-1}}{}{2} }.
\end{align*}

Noting  $F_{2,t}(w_t) = F_{1,t}(w_t) + v_t$ in the above relation, we have,

\begin{align}
\gamma_t \inprod{F_{2,t}(w_t)}{w_t-u} &+ \gamma_t\inprod{F_{1,t}(w_t)-F_{1,t}(u_{t-1})+ v_t}{u_t-w_t} \nonumber\\
&\quad\le \tfrac{1}{2}\bracket*{ \gnorm{u-u_{t-1}}{}{2} - \gnorm{u-u_{t}}{}{2} -\gnorm{u_t -w_t}{}{2}-\gnorm{w_t-u_{t-1}}{}{2} }.\label{eq:int_rel_NSEG_2}
\end{align}

Due to Lipschitz continuity of operator $F_{1,t}$ and Cauchy-Schwarz inequality, we have

\begin{align}
&\gamma_t\langle F_{1,t}(w_t)-F_{1,t}(u_{t-1})+ v_t,u_{t}-w_t\rangle
+\tfrac{1}{2}\bracket*{\gnorm{u_t -w_t}{}{2}+\gnorm{w_t-u_{t-1}}{}{2}} +\gamma_t^2\gnorm{v_t}{}{2} \nonumber\\
&\ge -\gamma_tL\gnorm{w_t-u_{t-1}}{}{}\gnorm{u_{t}-w_t}{}{}  -\gamma_t\gnorm{v_t}{}{}\gnorm{u_{t}-w_t}{}{} + \tfrac{1}{2}\gnorm{w_t-u_t}{}{2}+ \tfrac{1}{2}\gnorm{w_t-u_{t-1}}{}{2}+  \gamma_t^2\gnorm{v_t}{}{2} .\label{eq:int_rel_NSEG_1}
\end{align}

Note that we have

\begin{align*}
(\gamma_tL)^2 \le \tfrac{1}{2} \quad\Rightarrow\quad -\gamma_tL \gnorm{w_t-u_{t-1}}{}{}\gnorm{u_{t}-w_t}{}{}  + \tfrac{1}{4}\gnorm{w_t-u_t}{}{2} + \tfrac{1}{2}\gnorm{w_t-u_{t-1}}{}{2} \ge 0,\\
%\gamma_tL \le \tfrac{1}{4} \Rightarrow
-\gamma_t\gnorm{v_t}{}{}\gnorm{u_{t}-w_t}{}{} + \tfrac{1}{4}\gnorm{w_t-u_t}{}{2} + \gamma_t^2\gnorm{v_t}{}{2} \ge 0.
\end{align*}

Summing the above two relations along with (<ref>), we have

\[\gamma_t\langle F_{1,t}(w_t)-F_{1,t}(u_{t-1})+ v_t,u_{t}-w_t\rangle
+\tfrac{1}{2}\bracket*{\gnorm{u_t -w_t}{}{2}+\gnorm{w_t-u_{t-1}}{}{2}} +\gamma_t^2\gnorm{v_t}{}{2} \ge 0.\]

Using the above relation inside (<ref>),

we conclude the proof.

Now, we focus on the generalization properties of Algorithm <ref>. The next Theorem follows [22] to a
great extent. Its only major innovation is the stability analysis to control term (C) within the proof. This alternative analysis
shows that resampling the operator for the final update that defines $u_t$ is unnecesary to obtain optimal rates in the
nonprivate setting, which may be of independent interest.

Suppose all assumptions of Lemma <ref> are satisfied. Then, we have

\begin{eqnarray*}
\EE_{\bS,{\cal A}}\Big[\VIgap({\cal A}(\bS),F,\W)\Big] &\leq& \frac{1}{\Gamma_T}\Big( \gnorm{u-u_{0}}{}{2} +[17M^2+4(LD)^2+6MLD+2M]T\gamma^2 \\
\end{eqnarray*}

where $\Gamma_T = \tsum_{t=1}^T\gamma_t$, ${\cal A}(\bS)$ is the output of Algorithm <ref> on the dataset $\bS$ and expectation in the left hand  side is taken over the i.i.d. dataset draws, random sample batch choices, as well as noise $\bg^t_1, \bg^t_2$.
Moreover, if both noise schedule and stepsizes are constant throughout iterations, $\sigma_t\equiv\sigma$ and $\gamma_t\equiv\gamma$, then

\begin{eqnarray*}
\PP_{\bS,{\cal A}}\Big[\VIgap({\cal A}(\bS),F,\W) &>& \frac{1}{\Gamma_T}\Big( \gnorm{u-u_{0}}{}{2} +[17M^2+4(LD)^2+6MLD+2M]T\gamma^2 \\
&&+[8M+4M^2+4MLD]Td\sigma^2\gamma^2\Big) \Big]\leq 2\exp\{-dT/8\}.
\end{eqnarray*}

First note that by monotonicity of $F$,

\begin{equation}\label{eq:int_rel13}
\Gamma_T\langle F(u),\wb{w}_T-u \rangle
=\sum_{t=1}^{T} \gamma_t\langle F(u),w_t-u\rangle
\leq \sum_{t=1}^{T} \gamma_t\langle F(w_t),w_t-u\rangle.
\end{equation}

Let us define $\Delta_t := F(w_t) - F_{2,t}(w_t)$ and an auxiliary sequence
\begin{equation*}
z_{t}  := \proj_\W(z_{t-1}-\gamma_t\Delta_t)\mbox{, with } z_0=u_0, \label{eqn:def_prox_Delta}
\end{equation*}
for all $t \ge 1$. By definition of $z_t$, and the standard mirror-descent regret bound, we have

\begin{equation}\label{eq:int_rel14}
\sum_{t=1}^{T} \langle\gamma_t\Delta_t, z_t-u\rangle
\leq \frac12\gnorm{u-u_0}{}{2}-\frac12\gnorm{u-z_{T}}{}{2}+\sum_{t=1}^{T}\gamma_t^2\gnorm{\Delta_t}{}{2}.
\end{equation}

Noting that $\mathcal{A}(\bS) = \wb{w}_T$ is the output of Algorithm <ref> on dataset $\bS$ and using (<ref>),  we have,

\begin{align*}
\langle F(u),{\cal A}(\bS)-u\rangle
\leq \frac{1}{\Gamma_T}&\Big( \underbrace{\sum_{t=1}^{T}\gamma_t\langle F_{2,t}(w_t),w_t-u\rangle}_{A}+\underbrace{\sum_{t=1}^{T}\gamma_t\langle \Delta_t,z_t-u\rangle}_{B} +\underbrace{\sum_{t=1}^{T}\gamma_t\langle \Delta_t,w_t-z_t\rangle}_{C} \Big).%\\
%\leq\frac{1}{\Gamma_T}&\Big(\gnorm{u-u_0}{}{2}+\sum_{t=1}^{T}\gamma_t^2\Big(\big(4LD + M\gnorm{\bg_1^t}{}{}\big) \big[ LD+ M(1+\gnorm{\bg_2^t}{}{}) \big]+\gnorm{\Delta_t}{}{2}\Big)\\
%&+\sum_{t=1}^{T}\gamma_t \underbrace{\langle \Delta_t,w_t-y^t\rangle}_{=\mu_t} \Big),
\end{align*}

We now bound each of the terms above separately.

* Term A can be bounded by Lemma <ref> and the fact that $v_t = M(\boldsymbol{\xi}_t^2-\boldsymbol{\xi}_t^1)$.

* Term B can be bounded as in (<ref>).

* We proceed with analysis of  $\EE[C]$. Notice that neither $\|\Delta_t\|$ nor $\|w_t-z_t\|$ are guaranteed to be small (even in expectation). We will instead leverage the fact that our method is one-pass, together with an stability argument, to bound $C$.

Since $\langle  \Delta_t,w_t-z_t \rangle = \langle F(w_t)-F(w_t; \bB_t),w_t-z_t \rangle-M\langle  \boldsymbol{\xi}_t^2,w_t-z_t \rangle$, we will analyze both terms separately,

\begin{align*}
\EE\bracket{\langle \boldsymbol{\xi}_t^2,z_t -w_t\rangle} &= %\EE[\inprod{ \boldsymbol{\xi}_t^2}{u_{t-1}-w_t}] +
\EE[\inprod{ \boldsymbol{\xi}_t^2}{z_t-z_{t-1}}+\EE[\inprod{ \underbrace{\EE[\boldsymbol{\xi}_t^2|w_t, z_{t-1}]}_{=0} }{z_{t-1}-w_t} \\
&\leq \tfrac{M\gamma_t}{2}\|\boldsymbol{\xi}_t^2\|^2 +\tfrac{1}{2M\gamma_t}\|z_t-z_{t-1}\|^2 \\
&\leq \tfrac{M\gamma_t}{2}\|\boldsymbol{\xi}_t^2\|^2 +\tfrac{\gamma_t}{2M}\|\Delta_t\|^2 %\\
%&\leq& \frac{M\gamma_t}{2}\|\bg_2^t\|^2 +\frac{1}{M\gamma_t}[(\gamma_tM\|\bg_1^t\|)^2+\gamma_t^2\|M\bg_2^t-F_{\bB_t}(u_{t-1})\|^2]\\
%&\leq&  \frac{M\gamma_t}{2}\|\bg_2^t\|^2 +\gamma_t M[\|\bg_1^t\|^2+2\|\bg_2^t\|^2+2]\\
%&\leq& (3\|\bg_2^t\|^2+\|\bg_1^t\|^2+2)M\gamma_t.
\end{align*}

To bound the other term, we will use an stability argument. Let $\bB_{[1:t]}:=(\bB_1,\ldots,\bB_t)$ the sequence of batches used by the algorithm, and let $(\bB_{[1:t-1]},\bB_t^{\prime})=(\bB_1,\ldots,\bB_{t-1},\bB_t^{\prime})$. To make the dependence on the batches explicit, we also write $w_t=w_t(\bB_{[1:t]})$ (and analogously for $z_t)$. If
$\bB_T^{\prime}\sim{\cal P}^{B_t}$  is independent from $\bS$, then

\begin{align*}
\langle F&(w_t(\bB_{[1:t]}))-F(w_t(\bB_{[1:t]}) ; \bB_t),w_t(\bB_{[1:t]})-z_t(\bB_{[1:t]}) \rangle \\
&\leq \langle F(w_t(\bB_{[1:t]}))-F(w_t(\bB_{[1:t-1]},\bB_t^{\prime}); \bB_t),w_t(\bB_{[1:t]})-z_t(\bB_{[1:t]}) \rangle \\
&\quad+ \langle F(w_t(\bB_{[1:t-1]},\bB_t^{\prime}); \bB_t)-F(w_t(\bB_{[1:t]}); \bB_t),w_t(\bB_{[1:t]})-z_t(\bB_{[1:t]}) \rangle\\
&\leq \langle F(w_t(\bB_{[1:t]}))-F(w_t(\bB_{[1:t-1]},\bB_t^{\prime}); \bB_t),w_t(\bB_{[1:t]})-z_t(\bB_{[1:t]}) \rangle \\
&\quad + L\|w_t(\bB_{[1:t-1]},\bB_t^{\prime})-w_t(\bB_{[1:t]})\|\|w_t(\bB_{[1:t-1]})-z_t(\bB_{[1:t]})\| \\
&\leq  \langle F(w_t(\bB_{[1:t]}))-F(w_t(\bB_{[1:t-1]},\bB_t^{\prime}); \bB_t),w_t(\bB_{[1:t]})-z_t(\bB_{[1:t]}) \rangle  + 2MLD\gamma_t,
%&\leq&  \langle F_{\cal D}(w_t(\bB_{[1:t-1]},\bB_t^{\prime\prime}))-F_{\bB_t}(w_t(\bB_{[1:t-1]},\bB_t^{\prime})),w_t(\bB_{[1:t]})-z_t(\bB_{[1:t]}) \rangle \\
%&&+ L\|w_t(\bB_{[1:t-1]},\bB_t^{\prime}))-w_t(\bB_{[1:t]},\bB_t)\|\|w_t(\bB_{[1:t]})-z_t(\bB_{[1:t]})\| \\
\end{align*}

where in the last step we used that

\begin{align}
\|w_t&(\bB_{[1:t-1]},\bB_t^{\prime})-w_t(\bB_{[1:t]})\|\nonumber \\
&\leq \|u_{t-1}(\bB_{[1:t-1]})-\gamma_t F(u_t(\bB_{[1:t-1]}); \bB_t')
-u_{t-1}(\bB_{[1:t-1]})+\gamma_t F(u_t(\bB_{[1:t-1]}); \bB_t)\| \nonumber\\
&\leq 2M\gamma_t. \label{eqn:bd_stability_w}
\end{align}

Notice that an analog bound holds for $z_t$ as well.

On the other hand, if $\bB_t^{\prime\prime},\bB_t^{\prime\prime\prime}\sim{\cal P}^{B_t}$ are independent from $\bS,\bB_t^{\prime}$

\begin{align*}
&\langle F(w_t(\bB_{[1:t]}))-F(w_t(\bB_{[1:t-1]},\bB_t^{\prime}); \bB_t),w_t(\bB_{[1:t]})-z_t(\bB_{[1:t]}) \rangle \\
&= \langle F(w_t(\bB_{[1:t]}))-F(w_t(\bB_{[1:t-1]},\bB_t^{\prime}); \bB_t),w_t(\bB_{[1:t-1]}, \bB_t^{\prime\prime})-z_t(\bB_{[1:t-1]},\bB_t^{\prime\prime}) \rangle \\
& +2M[\|w_t(\bB_{[1:t-1]}, \bB_t^{\prime\prime})-w_t(\bB_{[1:t-1]}, \bB_t)\|+\|z_t(\bB_{[1:t-1]}, \bB_t^{\prime\prime})-z_t(\bB_{[1:t-1]}, \bB_t)\|]\\
&\leq \langle F(w_t(\bB_{[1:t]}))-F(w_t(\bB_{[1:t-1]},\bB_t^{\prime}); \bB_t),w_t(\bB_{[1:t-1]}, \bB_t^{\prime\prime})-z_t(\bB_{[1:t-1]},\bB_t^{\prime\prime}) \rangle +8M^2\gamma_t\\
&\leq \langle F(w_t(\bB_{[1:t-1]},\bB_t^{\prime\prime\prime}))-F(w_t(\bB_{[1:t-1]},\bB_t^{\prime}); \bB_t),w_t(\bB_{[1:t-1]}, \bB_t^{\prime\prime})-z_t(\bB_{[1:t-1]},\bB_t^{\prime\prime}) \rangle\\
&+  \langle F(w_t(\bB_{[1:t]}))- F(w_t(\bB_{[1:t-1]},\bB_t^{\prime\prime\prime})),w_t(\bB_{[1:t-1]}, \bB_t^{\prime\prime})-z_t(\bB_{[1:t-1]},\bB_t^{\prime\prime}) \rangle+8M^2\gamma_t\\
&\leq \langle F(w_t(\bB_{[1:t-1]},\bB_t^{\prime\prime\prime}))-F(w_t(\bB_{[1:t-1]},\bB_t^{\prime}); \bB_t),w_t(\bB_{[1:t-1]}, \bB_t^{\prime\prime})-z_t(\bB_{[1:t-1]},\bB_t^{\prime\prime}) \rangle\\
&+ L\|w_t(\bB_{[1:t]})-w_t(\bB_{[1:t-1]},\bB_t^{\prime\prime\prime})\| D+8M^2\gamma_t\\
&\leq \langle F(w_t(\bB_{[1:t-1]},\bB_t^{\prime\prime\prime}))-F(w_t(\bB_{[1:t-1]},\bB_t^{\prime}); \bB_t),w_t(\bB_{[1:t-1]}, \bB_t^{\prime\prime})-z_t(\bB_{[1:t-1]},\bB_t^{\prime\prime}) \rangle\\
&+ 2ML D\gamma_t+8M^2\gamma_t,
\end{align*}

where in the last step we used (<ref>). In particular, taking expectations
w.r.t. $\bB_t,\bB_t^{\prime},\bB_t^{\prime\prime},\bB_t^{\prime\prime\prime}$,

\begin{equation} \label{eqn:bd_martingale}
\EE_{\bB_t}[ \langle F(w_t(\bB_{[1:t]}))-F_{\bB_t}(w_t(\bB_{[1:t]})),w_t(\bB_{[1:t]})-z_t(\bB_{[1:t]}) \rangle|\bB_{[1:t-1]}] \leq  (2MLD + 8M^2)\gamma_t.
%|\langle F(w_t(\bB_{[1:t]}))-F_{\bB_t}(w_t(\bB_{[1:t]})),w_t(\bB_{[1:t]})-z_t(\bB_{[1:t]}) \rangle| &\leq&
\end{equation}

We combine the bounds above, to bound the SVI-gap

\begin{align}
\EE\bracket{\VIgap&({\cal A}(\bS),F,\W)} = \EE_{{\cal A},\bS}[\sup_{u \in \W} \inprod{F(u)}{\A(\bS) - u}]\nonumber\\
&\leq \EE\sup_{u \in \W}\Big[\tfrac{1}{\Gamma_T}\Big( \tsum_{t=1}^{T}\gamma_t\langle F_{2,t}(w_t),w_t-u\rangle+\tsum_{t=1}^{T}\gamma_t\langle \Delta_t,z_t-u\rangle + \tsum_{t=1}^{T}\gamma_t\langle \Delta_t,w_t-z_t\rangle \Big)\Big] \nonumber \\
&\leq \EE\Big[\tfrac{1}{\Gamma_T}\Big( \sup_{u \in \W}\big[\gnorm{u-u_{0}}{}{2} - \tfrac{1}{2}\gnorm{u-u_T}{}{2}
- \tfrac{1}{2}\gnorm{u-z_T}{}{2}\big] +\tsum_{t=1}^{T}\gamma_t^2\|\Delta_t\|^2+ \tsum_{t=1}^TM^2\gamma_t^2\gnorm{\boldsymbol{\xi}_t^2-\boldsymbol{\xi}_t^1}{}{2}\nonumber \\
%&&+ \sum_{t=1}^T[4(LD)^2+3MLD+(MLD+M^2/2)\|\bg_2^t\|^2+(MLD+M^2)(\|\bg_1^t\|^2+1)/2]\gamma_t^2 \nonumber \\
&\qquad\qquad + \tsum_{t=1}^{T}\gamma_t\langle \Delta_t,w_t-z_t\rangle \Big)\Big]\nonumber\\
&\leq\EE\big[ \tfrac{1}{\Gamma_T}\big( D^{2}+\tfrac{3}{2}\tsum_{t=1}^{T}2\gamma_t^2[(2M)^2+M^2\|\boldsymbol{\xi}_t^2\|^2]
%&+ \sum_{t=1}^T[4(LD)^2+3MLD+(MLD+ M^2/2)\|\bg_2^t\|^2+(MLD+M^2)(\|\bg_1^t\|^2+1)/2]\gamma_t^2 \nonumber\\
&\qquad\qquad+\tsum_{t=1}^TM^2\gamma_t^2\gnorm{\boldsymbol{\xi}_t^2-\boldsymbol{\xi}_t^1}{}{2} +\tsum_{t = 1}^TM^2\gamma_t^2\gnorm{\boldsymbol{\xi}_t^2}{}{2} \big)\big] \nonumber\\
&\leq \tfrac{1}{\Gamma_T}\big( D^2 +\tsum_{t=1}^{T}[20M^2+2MLD]\gamma_t^2 + \tsum_{t=1}^T5M^2\gamma_t^2\EE[\gnorm{\boldsymbol{\xi}_t^2}{}{2}]+\tsum_{t = 1}^TM^2\gamma_t^2\EE[\gnorm{\boldsymbol{\xi}_t^1}{}{2}]\big)%\nonumber\\
%&+\tsum_{t=1}^{T}(3M+MLD+ M^2)\|\bg_2^t\|^2\gamma_t^2+\sum_{t=1}^T(M+M^2+MLD)\|\bg_1^t\|^2\gamma_t^2\Big)
. \label{eqn:bound_mgf}
%&\leq& \frac{1}{\Gamma_T} \Big(D^2+\sum_{t=1}^T [2M+16M^2+2MLD]\gamma_t^2+\sum_{t=1}^T[6+2(LD)^2+5MLD+M^2]\sigma_t^2\gamma_t^2 \Big). \nonumber
\end{align}

Taking expectation, we obtain the first part of the result

\begin{align*}
\EE_{\bS}\big[\VIgap({\cal A}(\bS),F,\W)\big] &\leq \tfrac{1}{\Gamma_T}\big( D^{2} +[20M^2+2MLD]T\gamma^2+6M^2Td\sigma^2\gamma^2\big).
\end{align*}

Finally, for the high probability guarantee, first notice that if $\bg\sim{\cal N}(0,\sigma^2 I_{d\times d})$, then $\EE[\exp\{\lambda \bg_j^2/\sigma^2\}]=1/\sqrt{1-2\lambda}\leq\exp\{2\lambda^2+\lambda\}$ for $0<\lambda<1/4$, hence

\begin{eqnarray*}
\PP\Big(\sum_{t=1}^T\|\bg_1^t\|^2>2dT\sigma^2\Big) &= &\PP\Big(\sum_{t=1}^T\sum_{j\in[d]}\frac{\bg_{1,j}^2\lambda}{\sigma^2}>2dT\lambda\Big)
\,\,\leq \,\, \Big( \inf_{0<\lambda<1/4}\exp\{2\lambda^2-\lambda \} \Big)^{dT}
\,\, \leq \,\,  \exp\Big\{-\frac{dT}{8}\Big\}.
\end{eqnarray*}

The computation for $(\bg_2^t)_t$ is clearly analogous. These concentration bounds, combined
with (<ref>) provide the claimed high probability guarante,

\begin{align*}
\PP\big[\VIgap({\cal A}(\bS),F,\W) &>\tfrac{1}{\Gamma_T}\big( D^{2} +[20M^2+2MLD]T\gamma^2+12M^2Td\sigma^2\gamma^2\big) \big]\leq 2\exp\{-dT/8\}.
\end{align*}

Hence, we conclude the proof.

§.§ Differential privacy analysis of   method

We now proceed to establish the privacy guarantees for the single-pass variant of Algorithm <ref>. This is
a direct consequence of Propositions <ref> and <ref>, and the fact that
each operator evaluation has sensitivity bounded by $2M/B_t$.

Algorithm <ref> with batch sizes $(B_t)_{t\in[T]}$ and variance $\sigma_t^2=\frac{8M^2}{B_{t}^2}\frac{\ln(1/\privb)}{\priva^2}$
is $(\priva, \privb)$-differentially private.

We now apply the previous results to obtain population risk bounds for DP-SVI by the   method.

Algorithm <ref> with disjoint batches of size $B_t=B=\min\{\sqrt{d(\ln(1/\privb)}/\priva, n\}$, constant stepsize $\gamma_t\equiv\gamma=D/\big[M\sqrt{7T\big(1+\frac{8d}{B^2}\frac{\ln(1/\privb)}{\priva^2}\big)}\big]$ and variance
is $(\priva, \privb)$-differentially private and achieves $\VIgap{{\cal A}}{F}$ (for SVI) or $\SPgap{\A}{f}$ (for SSP) of

$$ O\Big(MD \max\Big\{ \frac{[d \ln(1/\privb)]^{1/4}}{\sqrt{n\priva}}, \frac{\sqrt{d\ln(1/\privb)}}{n\priva} \Big\}\Big). $$

Notice that in the corollary above, the gap is nontrivial iff $\sqrt{d\ln(1/\privb)}/[n\priva]<1$, which means that
the left hand side attains the max on the range where the gap is nontrivial.

Consider a SVI or SSP problem.
Let us recall that by Theorem <ref>, Algorithm <ref> achieves expected gap

$$ \frac{D^2}{\gamma T}+7M^2\gamma \Big(1+\frac{8d}{B^2}\frac{\ln(1/\privb)}{\priva^2} \Big).$$

Choosing $\gamma=D/\big[M\sqrt{7T\big(1+\frac{8d}{B^2}\frac{\ln(1/\privb)}{\priva^2}\big)}\big]$, we obtain

an expected gap

$$ \frac{2\sqrt{7}MD}{\sqrt{T}}\Big(1+\frac{\sqrt{8d}}{B\priva}\sqrt{\ln(1/\privb)}\Big)

where we used that for a single-pass algorithm, $n=2TB$ (this choice of $T$ exhausts the data when disjoint batches are chosen).

Recalling that $B=\min\{\sqrt{d\ln(1/\privb)}/\priva, n\}$. Then the expected gap is bounded by

$$ O\Big(MD \max\Big\{ \frac{[d\ln(1/\privb)]^{1/4}}{\sqrt{n\priva}}, \frac{\sqrt{d\ln(1/\privb)}}{n\priva} \Big\}\Big). $$

We observe that excess risk bounds of the same order for DP-SCO based on noisy SGD and the uniform stability of differential privacy have
been established [5].

Improving these bounds in DP-SCO required substantial efforts, which was only achieved recently
[4, 13, 3]. Furthermore, to the best of our knowledge, the upper bounds on the risk above are the first of their type for DP-SVI and DP-SSP, respectively. To improve upon them, we will follow the approach of [3], based on a multi-pass
empirical error convergence, combined with weak gap generalization bounds based on uniform stability.

Before stating the convergence results for Algorithm <ref>, we require a technical result stated below. This result was shown in [29] and we state and prove it here for the sake of completeness.

Let $F$ be any monotone Lipschitz operator. Consider the updates:

\begin{equation}\label{eqn:DPEG_iter1_new}
\begin{split}
x^t_{k-1} &\gets  \T{k-1}{\gamma_t  F} {x^{t-1}}, \\
x^t&\gets  \T{k}{\gamma_t  F} {x^{t-1}}.
\end{split}
\end{equation}

Consider the extragradient method in (<ref>) where $k = 2$. Suppose $\gamma_tL \le 1$. Then, for every $t=0,\ldots,T-1$, we have

$$ \langle \gamma_t F(x^t_{k-1}),x^t_{k-1}-u \rangle
\leq \frac12\|u-x^{t-1}\|^2-\frac12\|u-x^t\|^2%+2\gamma_t^2\gnorm{v_t}{}{2}
\quad(\forall u\in\W).$$

In what follows, let $D(u,v)=\frac12\|u-v\|^2$.
By the optimality conditions defining $x^t$, we have for all $u\in\W$

\begin{eqnarray}
&\langle \gamma_tF (x^t_{k-1}),x^t-u\rangle
&\leq D(u,x^{t-1})-D(u,x^t)-D(x^t,x^{t-1}) \notag\\
\Longrightarrow & \langle \gamma_tF(x^t_{k-1}),x^t_{k-1}-u\rangle
&\leq D(u,x^{t-1})-D(u,x^t)-D(x^t,x^{t-1}). \notag\\
\end{eqnarray}

Next, by the optimality conditions that define $x^t_{k-1}$,
we have

\begin{eqnarray*}
%\langle \gamma_t F(u_1^t),x^t_{k-1}-u_1^t\rangle &\leq&
%D(u_1^t,x^{t-1})-D(u_1^t,x^t_{k-1})-D(x^t_{k-1},x^{t-1}) \\
\langle \gamma_t F(x^{t-1}),x^t_{k-1}-x^t\rangle &\leq&
\end{eqnarray*}

Summing up the above bound with (<ref>),
we obtain

\begin{align*}
\langle \gamma_tF(x^t_{k-1}),x^t_{k-1}-u\rangle
%+\langle \gamma_t F(u_1^t),x^t_{k-1}-u_1^t\rangle
+\langle \gamma_t F(x^{t-1}),x_{k-1}^t-x^t\rangle \\
&\leq D(u,x^{t-1})-D(u,x^t)-D(x^t_{k-1},x^{t-1})-D(x^t,x_{k-1}^t).
\end{align*}

We note that, in order to prove the result, we need to show that

\begin{align}
\langle\gamma_tF(x^t_{k-1}),&x^t-x^t_{k-1}\rangle
+\langle \gamma_t F(x^{t-1}),x_{k-1}^t-x^t\rangle +D(x^t_{k-1},x^{t-1})+D(x^t,x_{k-1}^t)\geq 0.\label{eq:int_rel1}
\end{align}

Note that

\begin{align*}
\langle\gamma_tF(x^t_{k-1}),x^t-x^t_{k-1}\rangle
+\langle \gamma_t F(x^{t-1}),x_{k-1}^t-x^t\rangle &= \gamma_t\inprod{F(x^t_{k-1})- F(x^{t-1})}{x^t-x^t_{k-1}}\\
&\ge -\gamma_tL\gnorm{x^t_{k-1}-x^{t-1}}{}{}\gnorm{x^t-x^t_{k-1}}{}{}.
\end{align*}

Using (<ref>) along with above relation, we note that it suffices to show

\begin{align*}
&-\gamma_tL\gnorm{x^t_{k-1}-x^{t-1}}{}{}\gnorm{x^t-x^t_{k-1}}{}{} + \tfrac{1}{2}\gnorm{x^t_{k-1}-x^t}{}{2}+ \tfrac{1}{2}\gnorm{x^t_{k-1}-x^{t-1}}{}{2}\ge 0,
\end{align*}

which holds true in view of Young's inequality and the fact that $\gamma_t L \le 1$.
Hence, we conclude the proof.

Coming back to the proof of convergence for Algorithm <ref>, we show that the auxiliary sequence in (<ref>) satisfies a property that is crucial for convergence. Moreover, in view of Proposition <ref>, the same property is satisfied by the iterates of the algorithm, up to an additive error, since they are not too far from the auxiliary sequence.

We consider a variant of noisy extragradient method in (<ref>), (<ref>). Define:
\begin{align*}
F_{1,t}(\cdot) &:= F(\cdot; \xi_{t+1}) + M\bg_1,\\
F_{2,t}(\cdot) &:= F(\cdot; \xi_{t+1}) + M\bg_2.
\end{align*}

Consider the following algorithm:

\begin{eqnarray}
w_t &\gets&  T_{k-1}(u_t;\gamma_t  F_{1,t} ) \label{eqn:DPEG_iter1_new}\\
u_{t+1}&\gets & \prox_g(u_t-\gamma_t F_{2,t}(w_t)). \label{eqn:DPEG_iter2_new}
\end{eqnarray}

§ STABILITY OF NSEG AND OPTIMAL RISK FOR DP-SVI AND DP-SSP

The bounds established for DP-SVI are potentially suboptimal, and many of the past approaches used to
attain optimal rates for DP-SCO, such as privacy amplification by iteration, phased regularization, etc. appear to
encounter substantial barriers for their application to DP-SVI. In order to resolve this gap, we show
that for both DP-SVI and DP-SSP we can indeed obtain optimal rates,
which match those of DP-SCO. In order to achieve this we develop a multi-pass variant of the   method, which enjoys generalization performance due to its stability.

§.§ Stability of NSEG method

To analyze the stability of   it is useful to interpret the extragradient method as an approximation of the proximal point
algorithm. This connection has been established at least since [29]. Given a monotone
and $1$-Lipschitz operator $G:\bE\mapsto\bE$,

we define the $s$-extragradient
operator inductively as follows. First, $\T{0}{G}{\cdot}:\bE\mapsto\bE$ is defined as $\T{0}{G}{u}= \proj_{\W}(u)$. Then, for
$s\geq 0$

\begin{equation} \label{eqn:extrapolation}
\T{s+1}{G}{u}= \proj_{\W}(u- G(\T{s}{G}{u})).
\end{equation}

Given such operator, the (deterministic) extragradient method [24]
corresponds to,
starting from $u_0\in \W$, iterating

$$ u_{t+1} = \T{2}{\gamma F}{u_t} \qquad(\forall t\in[T-1]).$$

It is known that if $G$ is contractive, the recursion (<ref>) leads to
a fixed point $\T{}{G}{u}$, satisfying

\begin{equation} \label{eqn:fp_extrapolation}
\T{}{G}{u}= \proj_{\W}(u- G(\T{}{G}{u})).
\end{equation}

It is also easy to see that $\T{}{G}{\cdot}:\bE\mapsto\W$ is nonexpansive.

Let $F\in {\cal M}_{\W}^1(L,M)$ and $\W\subseteq \bE$ compact convex set with diameter $D>0$. Then, for all $s$
nonnegative integer, and $u,v\in\bE$,

\begin{equation} \label{eqn:fixed_pt_conv}
\| \T{s}{\gamma F}{u}- \T{}{\gamma F}{u}\| \leq
(\gamma L)^s\|\T{0}{\gamma F}{u}-\T{}{\gamma F}{u} \|
\end{equation}

\begin{equation} \label{eqn:near_nonexp_fixed_pt}
\| \T{s}{\gamma F}{u}- \T{s}{\gamma F}{v}\| \leq \|u-v\|+2D(\gamma L)^s.
\end{equation}

The first part, eqn. (<ref>), is proved by induction on $s$. The result clearly holds for $s=0$, and if $s\geq 1$, we use (<ref>) and (<ref>) to obtain

\begin{eqnarray*}
\| \T{s}{\gamma F}{u} - \T{}{\gamma F}{u} \| &=& \|\proj_{\W} (u-\gamma F(\T{s-1}{\gamma F}{u})- \proj_{\W}(u-\gamma F(\T{}{\gamma F}{u}))\| \\
&\leq&  \gamma\| F(\T{s-1}{\gamma F}{u})-F( \T{}{\gamma F}{u} )\| \,\,\leq\,\,  \gamma L\| \T{s-1}{\gamma F}{u}- \T{}{\gamma F}{u} \| \\
&\leq& (\gamma L)^s\| \T{0}{\gamma F}{u}-\T{}{\gamma F}{u}\|,
\end{eqnarray*}

where in the first inequality we used the nonexpansiveness of
the projection operator, next we used the $L$-Lipschitzness of $F$, and finally we used the inductive hypothesis to conclude.

The second part, eqn. (<ref>), is a direct consequence of (<ref>), the triangle inequality, and that $\T{0}{\gamma F}{u}$, $\T{}{\gamma F}{u}$, $\T{0}{\gamma F}{v}$, $\T{}{\gamma F}{v}\in \W$.

The next lemma shows an expansion upper bound for extragradient iterations. This type
of bound will be later used to establish the uniform argument stability of the NSEG algorithm.

Let $F_1,F_2:\bE\mapsto\bE$ monotone $L$-Lipschitz operators,
and $0\leq \gamma< 1/L$.
Let $u,v\in \W$, and $w,z,u^{\prime},v^{\prime}\in \W$ such that

\begin{eqnarray*}
w &=& \proj_{\W}(u-\gamma F_1(u)) %\T{k-1}{\gamma F_1}{u}
\qquad\qquad\,
z \,\,\,=\,\,  \proj_{\W}(v-\gamma F_1(v))  \\%\T{k-1}{\gamma F_1}{v}\\
u^{\prime} &=&  \proj_{\W}(u-\gamma F_2(w)) % \prox_{\gamma g}(u-\gamma F_2(w))
\qquad\qquad v^{\prime} \,\,=\,\,  \proj_{\W}(v-\gamma F_2(z)) .%\prox_{\gamma g}(v-\gamma F_2(z)).
\end{eqnarray*}

\begin{eqnarray}
%\mbox{ If $F_1=F_2$ } & & \|u^{\prime}-v^{\prime}\| \leq\|u-v\|+2L^2D \gamma^2  \label{eqn:exp_bd_F1eqF2}\\
%\mbox{ If $F_1\neq F_2$ } & &
\|w-z\| &\leq& \|u-v\|+2 L D\gamma,  \label{eqn:exp_bd_delta}\\
\|u^{\prime}-v^{\prime}\| &\leq& \|u-v\|+(\M_1+\M_2+2LD)L \gamma^2, \label{eqn:exp_bd_F1neqF2}
\end{eqnarray}

where $\M_1\triangleq \|F_1(\T{}{\gamma F_1}{u})-F_2(\T{}{\gamma F_2}{u})\|$ and $\M_2\triangleq \|F_1(\T{}{\gamma F_1}{v})-F_2(\T{}{\gamma F_2}{v})\|$.

By definition of $w$ and $z$, we have,

\begin{eqnarray*}
\|w-z\| &= & \| \proj_{\W}(u-\gamma F_1(u)) -\proj_{\W}(v-\gamma F_1(v))  \| \\
&\leq &\| \T{}{\gamma F_1}{u}- \T{}{\gamma F_1}{v}\|+\| \proj_{\W}(u-\gamma F_1(u))  - \T{}{\gamma F_1}{u}\|\\
&&+\| \proj_{\W}(v-\gamma F_1(v))  - \T{}{\gamma F_1}{v} \| \\
&\leq&\|u-v\|+ \gamma L [\| \T{0}{\gamma F_1}{u} - \T{}{\gamma F_1}{u} \|+\| \T{0}{\gamma F_1}{v}-T{}{\gamma F_1}{v}\|],
%2\gamma \|F(T_{k-2}^{\gamma F}(u))-F(T_{k-2}^{\gamma F}(v))\| \\
%&\leq& \|u-v\|+2(\gamma L)\|T_{k-2}^{\gamma F}(u)-T^{\gamma F}(v)\| \\
%&\leq&\|u-v\|+2 (\gamma L)^{k-1} \|T_{0}^{\gamma F}(u)-T_{0}^{\gamma F}(v)\|,
\end{eqnarray*}

where we used the nonexpansiveness of the operator
$\T{}{\gamma F_1}{\cdot }$ and Proposition <ref>. Moreover, since $u,v,\T{0}{\gamma F_1}{u}, \T{0}{\gamma F_1}{v} \in\W$, we have $\|w-z\|\leq \|u-v\|+2LD\gamma$, proving (<ref>).

Next, to prove (<ref>), we proceed as follows:

\begin{eqnarray*}
\|u^{\prime}-v^{\prime}\|
&= & \|\proj_{\W}(u-\gamma F_2(w))-\proj_{\W}(v-\gamma F_2(z))\|\\
&\leq& \| \T{}{\gamma F_2}{u} - \T{}{\gamma F_2}{v}\|+
\|\proj_{\W}(u-\gamma F_2(w))-\proj_{\W}(u-\gamma F_2(\T{}{\gamma F_2}{u}))\| \\
&& + \|\proj_{\W}(v-\gamma F_2(z))-\proj_{\W}(v-\gamma F_2(\T{}{\gamma F_2}{v}))\|\\
&\leq& \|u-v\|+\gamma L\|w- \T{}{\gamma F_2}{u}\|
+\gamma L \|z- \T{}{\gamma F_2}{v}\|.
\end{eqnarray*}

Using again Proposition <ref>, we have that

\begin{eqnarray*}
\|w- \T{}{\gamma F_2}{u} \|
&=& \| \T{1}{\gamma F_1}{u}- \T{}{\gamma F_2}{u} \| \\ &\leq& \| \T{1}{\gamma F_1}{u}-\T{}{\gamma F_1}{u} \| + \| \T{}{\gamma F_1}{u} - \T{}{\gamma F_2}{u} \|  \\
&\leq&  L D\gamma+\|\proj_{\W}\big(u-\gamma F_1( \T{}{\gamma F_1}{u})\big)-\proj_{\W}\big(u-\gamma F_2( \T{}{\gamma F_2}{u})\big)\| \\
&\leq&  L D\gamma+\gamma \|F_1(\T{}{\gamma F_1}{u})-F_2(\T{}{\gamma F_2}{u})\|\\
&\leq&  L D\gamma+\M_1\gamma.
\end{eqnarray*}

An analog bound can be obtained for $\|z-\T{}{\gamma F_2}{v}\|$:

$$ \|z-\T{}{\gamma F_2}{v}\|\leq LD\gamma+\M_2\gamma, $$

concluding the claimed bound (<ref>):

$$ \|u^{\prime}-v^{\prime}\|\leq \|u-v\|+ L[\M_1+\M_2+ 2LD]\gamma^2.$$

The Expansion Lemma above allows us to bound how much would two trajectories of the   method may deviate, given two pairs of sequences of operators $F_{1,t},F_{2,t}$ and $F_{1,t}^{\prime},F_{2,t}^{\prime}$.
The bounds we will obtain from this analysis will give us direct bounds on the UAS
for the NSEG  method.

Let $F_{1,t},F_{2,t}$ and $F_{1,t}^{\prime},F_{2,t}^{\prime}$ be $L$-Lipschitz operators,

and $0\leq \gamma_t< 1/L$ for all $t\in[T]$. Let
$\{(u_t,w_t)\}_{t\in[T]}$ and $\{(v_t,z_t)\}_{t\in[T]}$ be the sequences resulting from Algorithm <ref>, with operators $\{(F_{1,t},F_{2,t})\}_{t\in[T]}$ and $\{(F_{1,t}^{\prime},F_{2,t}^{\prime})\}_{t\in[T]}$, respectively; and starting from $u^0=v^0$. Let

\begin{eqnarray*}
\Delta_{1,t} &\triangleq& \sup_{u\in\W}\|F_{1,t}(u)-F_{1,t}^{\prime}(u)\|, \\
\Delta_{2,t}  &\triangleq&\sup_{u\in\W}\|F_{2,t}(u)-F_{2,t}^{\prime}(u)\|,\\
\M_{1,t}       &\triangleq& \|F_{1,t}(\T{}{\gamma F_{1,t}}{u_{t-1}})-F_{2,t}(\T{}{\gamma F_{2,t}}{u_{t-1}})\|, \mbox{ and } \\
\M_{2,t}       &\triangleq& \|F_{1,t}(\T{}{\gamma F_{1,t}}{v_{t-1}})-F_{2,t}(\T{}{\gamma F_{2,t}}{v_{t-1}})\|;
\end{eqnarray*}

then, for al $t=0,\ldots,T$,

\begin{eqnarray}
\nu_t \triangleq \|u_t-v_t\| &\leq& \sum_{s=1}^{t}\Big( [\M_{1,t}+\M_{2,t}+2LD]L\gamma_s^2+L\Delta_{1,s}\gamma_s^2+\Delta_{2,s}\gamma_s\Big)
%\sum_{s=0}^{t-1} \Big( [4M+2LD]L\gamma_s^2+2(L\gamma_s)^2\Delta_{1,s}+2L\gamma_s\Delta_{2,s} \Big);
\label{eqn:bd_nu_F1neqF2}\\
\delta_t \triangleq  \|w_t-z_t\| &\leq&  \sum_{s=1}^{t-1}\Big( [\M_{1,t}+\M_{2,t}+2LD]L\gamma_s^2+L\Delta_{1,s}\gamma_s^2+\Delta_{2,s}\gamma_s\Big)
%\sum_{s=0}^{t-1}\Big( [4M+2LD]L\gamma_s^2+\Delta_{2,s}\gamma_s\Big)
%+ \sum_{s=0}^{t}L\Delta_{1,s}\gamma_s^2 +2LD\gamma_t
\label{eqn:bd_delta_F1neqF2}
\end{eqnarray}

Clearly, $\nu_0=0$. Let us now derive a recurrence for both $\nu_t$ and $\delta_t$.

\begin{eqnarray*}
\delta_{t} &=& \| \T{1}{\gamma_t F_{1,t}}{u_{t-1}} - \T{1}{\gamma_t F_{1,t}^{\prime}}{v_{t-1}} \|  \\
&\leq&  \| \T{1}{\gamma_t F_{1,t}}{u_{t-1}} - \T{1}{\gamma_t F_{1,t}}{v_{t-1}} \|+
\| \T{1}{\gamma_t F_{1,t}}{v_{t-1}} - \T{1}{\gamma_t F_{1,t}^{\prime}}{v_{t-1}} \| \\
&\leq& \|u_{t-1}-v_{t-1}\|+2LD\gamma_t + \| \T{1}{\gamma_t F_{1,t}}{v_{t-1}} - \T{1}{\gamma_t F_{1,t}^{\prime}}{v_{t-1}}\|,
\end{eqnarray*}

where in the last inequality we used inequality (<ref>).
Let us bound now the rightmost term above,

\begin{eqnarray}
\| \T{1}{\gamma_t F_{1,t}}{v_{t-1}} - \T{1}{\gamma_t F_{1,t}^{\prime}}{v_{t-1}} \|
&=& \|\proj_{\W}(v_{t-1}-\gamma_t F_{1,t}(v_{t-1}))-\proj_{\W}(v_{t-1}-\gamma_t F_{1,t}^{\prime}(v_{t-1}))\|  \notag\\
&\leq&\gamma_t \| F_{1,t}(v_{t-1})- F_{1,t}^{\prime}(v_{t-1})\| \notag \\
&\leq & \Delta_{1,t}\gamma_t.
\label{eqn:rec_bd_stab}
\end{eqnarray}

We conclude that

\begin{equation} \label{eqn:bd_delta}
\delta_{t} \leq \nu_{t-1} + 2LD\gamma_t + \Delta_{1,t}\gamma_t.
\end{equation}

\begin{align*}
\nu_{t} =& \|u_{t}-v_{t}\| \leq \|\proj_{\W}(u_{t-1}-\gamma_t F_{2,t}(w_t))-\proj_{\W}(v_{t-1}-\gamma_t F_{2,t}^{\prime}(z_t))\| \notag \\
=& \|\proj_{\W}(u_{t-1}-\gamma_t F_{2,t}(w_{t}))-\proj_{\W}(v_{t-1}-\gamma_t F_{2,t}(z_t))\|\notag \\
&+\|\proj_{\W}(v_{t-1}-\gamma_t F_{2,t}(z_t))-\proj_{\W}(v_{t-1}-\gamma_t F_{2,t}^{\prime}(z_t))\| \notag\\
\leq& \|\proj_{\W}(u_{t-1}-\gamma_t F_{2,t}( \T{1}{\gamma_t F_{1,t}}{u_{t-1}}))-\proj_{\W}(v_{t-1}-\gamma_t F_{2,t}(\T{1}{\gamma_t F_{1,t}^{\prime}}{v_{t-1}}))\| \notag\\
&+\gamma_t \|F_{2,t}(z_t)-F_{2,t}^{\prime}(z_t)\| \notag \\
\leq&\|\proj_{\W}(u_{t-1}-\gamma_t F_{2,t}(\T{1}{\gamma_t F_{1,t}}{u_{t-1}}))-\proj_{\W}(v_{t-1}-\gamma_t F_{2,t}(\T{1}{\gamma_t F_{1,t}}{v_{t-1}}))\| \notag \\
&+\|\proj_{\W}(v_{t-1}-\gamma_t F_{2,t}(\T{1}{\gamma_t F_{1,t}}{v_{t-1}}))-\proj_{\W}(v_{t-1}-\gamma_t F_{2,t}(\T{1}{\gamma_t F_{1,t}^{\prime}}{v_{t-1}}))\| +\gamma_t \Delta_{2,t}  \notag \\
\mleq{(i)}& \|u_{t-1}-v_{t-1}\|+[\M_{1,t}+\M_{2,t}+2LD]L\gamma_t^2
+\gamma_tL\| \T{1}{\gamma_t F_{1,t}}{v_{t-1}} - \T{1}{\gamma_t F_{1,t}^{\prime}}{v_{t-1}}\|+\gamma_t \Delta_{2,t} %\label{use:lem_exp}
\\
\mleq{(ii)}& \nu_{t-1} + [\M_{1,t}+\M_{2,t}+2LD]L\gamma_t^2+L\Delta_{1,t}\gamma_t^2+\Delta_{2,t}\gamma_t, %\label{use:rec_bd_stab}
\end{align*}

where in inequality (i),
we used Lemma <ref> (more precisely, inequality (<ref>)),
and in inequality (ii),
we used (<ref>). Unraveling the above recursion, we get that for all $t\in[T]$,

$$ \nu_t \leq \sum_{s=1}^{t}\Big( [\M_{1,t}+\M_{2,t}+2LD]L\gamma_s^2+L\Delta_{1,s}\gamma_s^2+\Delta_{2,s}\gamma_s\Big).$$

Finally, we combine the bound above with (<ref>), to conclude that for all $t\in[T]$:

$$ \delta_t \leq \sum_{s=1}^{t-1}\Big( [\M_{1,t}+\M_{2,t}+2LD]L\gamma_s^2+L\Delta_{1,s}\gamma_s^2+\Delta_{2,s}\gamma_s\Big)

The next theorem provides in-expectation and high probability upper bounds for the   method. Despite the fact that we will not particularly apply the latter bounds, we believe these may be of independent interest.

The   method (Algorithm <ref>) for closed and convex domain $\W\subseteq \bE$
with diameter $D$, operators in ${\cal M}_{\W}^1(L, M)$ and stepsizes $0<\gamma_t\leq 1/L$, satisfies the
following uniform argument stability bounds:

* Batch method $\AbatchEG$: Where given dataset $\bS$, $F_{1,t}=F_{\bS}+\bg_1^t$, and $F_{2,t}=F_{\bS}+\bg_2^t$. In expectation,

$$\sup_{\bS\simeq \bS^{\prime}}\EE_{\AbatchEG}[\delta_{\AbatchEG}(\bS,\bS^{\prime})]\leq
\sum_{t=0}^{T-1}\Big( [4M+2LD+4\sqrt{d}\sigma]L\gamma_t^2+\frac{2ML}{n}\gamma_t^2+\frac{2M}{n}\gamma_t\Big)
+ \frac{1}{T}\sum_{t=1}^{T}\big(\frac{2ML}{n}+2LD\big)\gamma_t,$$

And for constant stepsize $\gamma_t\equiv\gamma$,
there exists a universal constant $K>0$, such that for any $0<\theta<1$, with probability $1-\theta$:

\begin{equation} \label{eqn:h_p_bd_batch_NSEG}
\sup_{\bS\simeq \bS^{\prime}}\delta_{\AbatchEG}(\bS,\bS^{\prime})
\leq 4[T\sqrt{d}\sigma+\sigma\sqrt{Kd\ln(1/\theta)}]L\gamma^2+
+ \big(\frac{2ML}{n}+2LD\big)\gamma.
\end{equation}

* Sampled with replacement $\AreplEG$: Where given dataset $\bS$, $F_{1,t}=F_{\bbeta_{i(1,t)}}+\bg_1^t$, and $F_{2,t}=F_{\bbeta_{i(2,t)}}+\bg_2^t$, for $i(1,t),i(2,t)\sim\mbox{Unif}([n])$, independently. In expectation,

\begin{equation} \label{eqn:exp_UAS_stoch_EG}
\sup_{\bS\simeq \bS^{\prime}}\EE[\delta_{\AreplEG}(\bS,\bS^{\prime})]
\,\,\leq\,\, \sum_{t=0}^{T-1}\Big( [4M+2LD+4\sqrt{d}\sigma]L\gamma_t^2+\frac{2 ML}{n}\gamma_t^2+\frac{2 M}{n}\gamma_t\Big)
+ \frac{1}{T}\sum_{t=1}^{T}\big(\frac{2ML}{n}+2LD\big)\gamma_t.
\end{equation}

And for constant stepsize $\gamma_t\equiv\gamma$, there exists a universal constant $K>0$, such that
for any $0<\theta<1/[2n]$, with probability $1-\theta$:

\begin{multline}  \label{eqn:HP_UAS_stoch_EG}
\sup_{\bS\simeq \bS^{\prime}}\delta_{\AreplEG}(\bS,\bS^{\prime})\leq  4[T\sqrt{d}\sigma+\sigma\sqrt{Kd\ln(2/\theta)}]L\gamma^2+ [4M+2LD]LT\gamma^2 +2LD\gamma \\
\end{multline}

Let $\bS\simeq \bS^{\prime}$. Then

* Batch method. Notice that for the batch case $F_{1,t}=F_{\bS}+\bg_1^t$, and
$F_{1,t}^{\prime}=F_{\bS^{\prime}}+\bg_1^t$; and $F_{2,t}=F_{\bS}+\bg_2^t$, and
Then, it is easy to see that $\Delta_{1,t}\leq2M/n$ and $\Delta_{2,t}\leq2M/n$. On the other hand,
since the operators are $M$ bounded and since noise addition is Gaussian

\begin{eqnarray}
\EE[\M_{1,t}] &=& \EE[\|F_{1,t}(\T{}{\gamma F_{1,t}}{u_{t-1}})-F_{2,t}(\T{}{\gamma F_{2,t}}{u_{t-1}})\|] \nonumber \\
&\leq&  \EE[\|F_{\bB_t^1}(\T{}{\gamma F_{1,t}}{u_{t-1}})+\bg_1^t\|]+\EE[\|F_{\bB_t^2}(\T{}{\gamma F_{2,t}}{u_{t-1}})+\bg_2^t\|] \nonumber\\
&\leq& 2M+\EE[\|\bg_1^t\|+\|\bg_2^t\|] \leq 2[M+\sqrt{d}\sigma], \label{eqn:bd_noise_stab}
\end{eqnarray}

and an analog bound holds for $\EE[\M_{2,t}]$. Hence, by Lemma <ref>:

\sum_{t=0}^{T-1}\Big( [4M+2LD+4\sqrt{d}\sigma]L\gamma_t^2+\frac{2ML}{n}\gamma_t^2+\frac{2M}{n}\gamma_t\Big)
+ \frac{1}{T}\sum_{t=1}^{T}\big(\frac{2ML}{n}+2LD\big)\gamma_t,$$

which proves the claimed bound.

For the high probability bound, we use that the norm of a Gaussian vector is $Kd\sigma^2$-subgaussian,
for a universal constant $K>0$ (see, e.g. <cit.>), and therefore
$\EE[\exp\{\lambda(\|\bg_i^t\|-\sigma\sqrt{d})\}]\leq \exp\{Kd\sigma^2\lambda^2\}$; hence by the Chernoff-Crámer bound,
for any $\alpha>0$

\begin{eqnarray*}
\PP\Big[\sum_{t\in[T]} \big(\|\bg_1^t\|+\|\bg_2^t\|\big) >(2+\alpha)T\sqrt{d}\sigma)\Big]
&\leq& \exp\{-\lambda\alpha T\sqrt{d}\sigma\}\Big( \exp\{2Kd\sigma^2\lambda^2\} \Big)^T\\
&=& \exp\{T(2Kd\sigma^2\lambda^2-\alpha\sqrt{d}\sigma\lambda)\}.
\end{eqnarray*}

Choosing $\lambda=\alpha/[4K\sqrt{d}\sigma]$ and $\alpha=\frac{2\sqrt{K}}{T}\sqrt{\ln(1/\theta)}$, we get

\begin{equation} \label{eqn:concentration_gaussians}
\PP\Big[\sum_{t\in[T]} \big(\|\bg_1^t\|+\|\bg_2^t\|\big) > 2T\sqrt{d}\sigma+2\sigma\sqrt{Kd\ln(1/\theta)}\Big]\leq\theta.
\end{equation}

This guarantee, together with the rest of the terms appearing in our previous stability bound (which hold w.p. 1)
proves (<ref>).

* Sampled with replacement.

Let $i\in[n]$ be the coordinate where $\bS$ and $\bS^{\prime}$ may differ.
Let $i(1,t),i(2,t)\sim \mbox{Unif}([n])$ i.i.d., for $t\in [T]$. Now we apply Lemma <ref> with
$F_{1,t}=F_{\bbeta_{i(1,t)}}+\bg_1^t$, and $F_{1,t}^{\prime}=F_{\bbeta_{i(1,t)}^{\prime}}+\bg_1^t$;
and $F_{2,t}=F_{\bbeta_{i(2,t)}}+\bg_2^t$, and $F_{2,t}^{\prime}=F_{\bbeta_{i(2,t)}^{\prime}}+\bg_2^t$.
Hence we have that $(\Delta_{1,t})_{t\in[T]}$ and $(\Delta_{2,t})_{t\in[T]}$ are
sequences of independent r.v. with expectation bounded by $2M/n$. Therefore,
by Lemma <ref>  (eqn. (<ref>)), and following the
steps that lead to inequality (<ref>), we have:

\begin{eqnarray*}
\EE[\delta_{\cal A}(\bS,\bS^{\prime})]
&\leq& \sum_{t=0}^{T-1}\Big( [4M+2LD+4\sqrt{d}\sigma]L\gamma_t^2+\frac{2ML}{n}\gamma_t^2+\frac{2M}{n}\gamma_t\Big)
+ \frac{1}{T}\sum_{t=1}^{T}\big(\frac{2ML}{n}+2LD\big)\gamma_t.
%\sum_{t=1}^{T}\Big( [4M+2LD]L\gamma_t^2+\frac{2ML}{n}\gamma_t^2\Big)+\frac{2M}{n}\gamma_t\Big) +\frac{2LD}{T}\sum_{t=1}^T\gamma_t .
\end{eqnarray*}

Finally for the high-probability bound, note that for any realization of the algorithm randomness, we have

$$ \delta_{\cal A}(\bS,\bS^{\prime}) \leq \sum_{t=1}^T[4M+2(\|\bg_1^t\|+\|\bg_2^t\|)+2LD]L\gamma_t^2+\frac{2LD}{T}\sum_{t=1}^T\gamma_t

We additionally assume constant stepsize, $\gamma_t\equiv \gamma>0$.

Hence, we can resort on concentration of sums of Bernoulli random variables, which guarantees that

$$ \PP\Big[ \sum_{t=1}^T\Delta_{1,t} > (1+3\log(2/\theta))\frac{2MT}{n} \Big]
\leq \exp\{-\log(2/\theta)\}=\frac{\theta}{2}. $$

An analog bound can be established for $\Delta_{2,t}$, which together with bound (<ref>)
leads to

\begin{multline*}
\PP_{\AreplEG}\Big[ \delta_{\AreplEG}(\bS,\bS^{\prime})> 4[T\sqrt{d}\sigma+\sigma\sqrt{Kd\ln(1/\theta)}]L\gamma^2+ [4M+2LD]LT\gamma^2 +2LD\gamma \\
+\big(1+3\log\big(\frac{2}{\theta}\big)\big)\frac{2MT}{n}(L\gamma^2+\gamma/T+\gamma) \Big]\leq 2\theta.
\end{multline*}

Notice this bound only depends on our choice of $i$, and it is otherwise uniform over all $\bS\simeq \bS^{\prime}$.
Finally, by a union bound on $i\in[n]$ (together with a renormalization of $\theta$), we have that

\begin{multline*}
\PP_{\AreplEG}\Big[ \sup_{\bS\simeq \bS^{\prime}} \delta_{\AreplEG}(\bS,\bS^{\prime})> 4[T\sqrt{d}\sigma+\sigma\sqrt{Kd\ln(2/\theta)}]L\gamma^2+ [4M+2LD]LT\gamma^2 +2LD\gamma \\
+\big(1+3\log\big(\frac{4n}{\theta}\big)\big)\frac{2MT}{n}(L\gamma^2+\gamma/T+\gamma) \Big]\leq \theta.
\end{multline*}

§.§ Optimal risk for DP-SVI and DP-SSP by  method

Now we use our stability and risk bounds for NSEG to derive optimal risk bounds for DP-SSP. For this, we use
the sampled with replacement variant, $\AreplEG$.

\begin{equation} \label{eqn:NSEG_setting}
F_{1,t}(\cdot)=F_{\boldsymbol{\beta}_{i(1,t)}}(\cdot)+\bg_1^t; \qquad
%F_{1,t}(\cdot)= F(\cdot\,;S_{i(t)})+M\boldsymbol{\xi}_1^t; \qquad
%F_{2,t}(\cdot)= F(\cdot\,;S_{i(t)})+M\boldsymbol{\xi}_2^t
%\quad \Big(\,i(t)\sim\mbox{Unif}([n]),\,\,\boldsymbol{\xi}_1^t,\boldsymbol{\xi}_2^t\sim {\cal N}(0,\sigma_t^2 )\,\Big).
\end{equation}

Using the moments accountant method (Theorem <ref>) one can show the following.

Algorithm <ref> with operators given by eqn. (<ref>) and
$\sigma_t^2=8M^2\log(1/\privb)/\priva^2$, is $(\priva,\privb)$-differentially private.

Consider an instance of the  (<ref>) or (<ref>) problem.

Let $\A$ be the sampled with replacement variant (<ref>) of  method (Algorithm <ref>), with
$\gamma_t=\gamma=\min\{D/M,1/L\}/[n\max\{\sqrt n, \sqrt{d \ln(1/\privb)}/\priva\}]$, $\sigma_t^2=8M^2\log(1/\privb)/\priva^2$, $T=n^2$.
Then, $\wVIgap{\A}{F}$ (for SVI) or $\wSPgap{\A}{f}$ (for SSP) are bounded by

O\Big( (MD+LD^2)\max\Big\{\frac{1}{\sqrt n}, \frac{\sqrt{d\ln(1/\privb)}}{n\priva}\Big\} +\frac{MLD}{n^{5/2}}\Big).

Notice that assuming $n=\Omega(\min\{\sqrt{L},\sqrt{M/D}\})$ the bound of the Theorem simplifies to

O\Big(  (MD+LD^2)\max\Big\{\frac{1}{\sqrt n}, \frac{\sqrt{d\ln(1/\privb)}}{n\priva}\Big\} \Big).

This is quite a mild sample size requirement. In this range, when $M\geq LD$, our upper bound matches the excess risk bounds for DP-SCO [4],
and we will show these rates are indeed optimal for DP-SVI and DP-SSP as well

Given that our bounds for SVI and SSP are analogous, we proceed indistinctively for both problems.

First, let us bound the empirical accuracy of the method. By Theorem <ref>, together with the fact
that sampling with replacement is an unbiased stochastic oracle for the empirical operator:

\begin{eqnarray*}
\EE_{\bS}\Big[\mbox{EmpGap}(\A,\bS)\Big]
&\leq& \frac{1}{\gamma T}\Big( \frac{D^2}{2}+7M^2T\gamma^2(1+\frac{8d \log(1/\privb)}{\priva^2}) \Big)\\
&\leq& \frac{D^2}{2n}\max\{M/D,L\} \max\{\sqrt{n},\sqrt{d\ln(1/\privb)}/\varepsilon\}
+\frac{7M^2\min\{D/M,1/L\}}{n\max\{\sqrt{n},\sqrt{d\ln(1/\privb)}/\varepsilon\}} \frac{9d\ln(1/\privb)}{\varepsilon^2} \\
%&\leq& \frac{MD}{2}\max\Big\{\frac{1}{\sqrt n}, \frac{\sqrt{d\ln(1/\privb)}}{n\priva}\Big\}+7MD\frac{9\sqrt{d\ln(1/\privb)}}{n\priva}\\
&=& O\Big( (MD+LD^2)\max\Big\{\frac{1}{\sqrt n}, \frac{\sqrt{d\ln(1/\privb)}}{n\priva}\Big\} \Big),
\end{eqnarray*}

where $\mbox{EmpGap}(\A,\bS)$ is $\EmpVIgap{\A}{F_{\bS}}$ or
$\EmpSPgap{\A}{f_{\bS}}$ for an SVI or SSP problem, respectively.

Next, by Theorem <ref>, we have that ${\cal A}$ (or $x(\bS)$ and $y(\bS)$, for the SSP case)
are UAS with parameter

\begin{eqnarray*}
%\Delta_x+\Delta_y &\leq&
\delta&=& [4M+2LD+4\sqrt{d}\sigma]LT\gamma^2+\frac{2 ML}{n}T\gamma^2+\frac{2 M}{n}T\gamma
+ \big(\frac{2ML}{n}+2LD\big)\gamma\\
&\leq& \Big(\frac{4LD^2}{M}+2D\Big)\frac1n+\frac{8LD^2\sqrt{2d\ln(1/\eta)}}{M\varepsilon n}+\frac{2LD^2}{M}\frac{1}{n^{3/2}}+\frac{2D}{\sqrt n}+\frac{2LD}{n^{5/2}}+\frac{2D}{n^{3/2}}\\
&=& O\Big(\frac{1}{M}\cdot\Big(\frac{MD+LD^2}{n}+\frac{LD^2\sqrt{d\ln(1/\eta)}}{\varepsilon n} +\frac{MLD}{n^{5/2}}\Big) \Big).
%\Big( \frac{4LD^2/M+2L^2D^3/M^2}{n}+\frac{8\sqrt2LD^2}{M}\frac{\sqrt{d\ln(1/\privb)}}{\priva n}+\frac{2LD^2}{Mn^2}+\frac{2D}{\sqrt{n}}
%+2\big(\frac{LD}{n}+\frac{LD^2}{M}\big)\frac{1}{n^{3/2}} \Big).
\end{eqnarray*}

Hence, noting that empirical risk upper bounds weak empirical gap and using Proposition <ref> or Proposition <ref> (depending on whether
the problem is an SSP or SVI, respectively), we have that the risk is upper bounded by its empirical risk plus $M\delta$, where
$\delta$ is the UAS parameter of the algorithm; in particular,   is bounded by

\begin{eqnarray*}
\wVIgap{\A}{F} &\le& \EE_{\bS}[\EmpVIgap{\A}{F_{\bS}}]+M\delta\\
&=& O\Big( (MD+LD^2)\max\Big\{\frac{1}{\sqrt n}, \frac{\sqrt{d\ln(1/\privb)}}{n\priva}\Big\} +\frac{MLD}{n^{5/2}}\Big),
%O\Big(  (1+LD/M)\frac{LD^2}{n}+(1+L/n^2)\frac{MD}{\sqrt{n}}+(MD+LD^2)\frac{\sqrt{d\ln(1/\privb)}}{\priva n} \Big),
\end{eqnarray*}

Similar claims can be made $\wSPgap{\A}{f_\bS}$. Hence, we conclude the proof.

§ THE NOISY INEXACT STOCHASTIC PROXIMAL POINT METHOD

In the previous sections, we presented   method with its single-pass and multipass variants and provided optimal risk guarantees for DP-SVI and DP-SSP problems in $O(n^2)$ stochastic operator evaluations. In the rest of the paper, our aim is to provide another algorithm that can achieve the optimal risk for both of these problems with much less computational effort. Towards that end, consider the following algorithm:
Noisy Inexact Stochastic Proximal Point (NISPP) Method

[1]
Input: $w_0 \in \W$
$k = 0, 1, \dots, K$
Sample a batch $\bB_{k+1} \subseteq \bS$.
$u_{k+1} \gets \text{VI}_\nu(\W, F_{\bB_{k+1}}(\cdot)+ \lambda_k (\cdot-w_k))$.

$w_{k+1} \gets u_{k+1} + \boldsymbol{\xi}_{k+1}$, where $\boldsymbol{\xi}_{k+1} \sim \mathcal{N}(0, \sigma_{k+1}^2 \mathbb{I}_d)$

$\wb{w}_K :=  (\tsum_{k =0}^K \gamma_kw_{k+1})/(\tsum_{k = 0}^K\gamma_k)$
Output: $\proj_{\W}(\wb{w}_K)$

In the above algorithm, we leave few things unspecified which will be stated later during convergence and privacy analysis. Here, we detail some key features of the above algorithm. In line 3, we sample a batch $\bB_{k+1}$ of size $B_{k+1}=|\bB_{k+1}|$.

Similar to the , we will look at two different variants of NISPP method: single-pass and multi-pass. Depending on the type of the method, we will specify the sampling mechanism.
In line 4 of Algorithm <ref>, we have $u_{k+1}$ is a $\nu$-approximate strong VI solution of the mentioned VI problem for some $\nu \ge 0$, i.e.,

\begin{equation}
\label{eq:int_rel_prox_1}
\inprod{F_{\bB_{k+1}}(u_{k+1}) + \lambda_k (u_{k+1}-w_k)}{w-u_{k+1}} \ge -\nu \qquad (\forall w \in \W).
\end{equation}

Note that if $\nu = 0$ then this is an exact solution satisfying (<ref>) with operator $F$ replaced by $F(\cdot) + \lambda_k(\cdot- w_k)$. For $\nu>0$, we obtain that $u_{k+1}$ is an inexact solution satisfying solution criterion up to $\nu$ additive error. In line 5, we add a Gaussian noise to $u_{k+1}$ in order to preserve privacy. The resulting iterate $w_{k+1}$ can be potentially outside the set $\W$. Hence, in line 7, the ergodic average $\wb{w}_K$ can be outside $\W$. In order to preserve feasibility of the solution, we project $\wb{w}_K$ onto set $\W$ and output it as a solution in line 8. Projection of the average in line 8, as opposed to projection individual $w_{k+1}$ in line 5 is crucial for convergence guarantee of Algorithm <ref>.

In the rest of this section, we exclusively deal with the single-pass version of NISPP method,

i.e., we assume that batches $\{\bB_{k+1}\}_{k=0,\ldots,K}$ are disjoint subsets of the dataset $\bS$.

We start with the convergence guarantees of single-pass NISPP method.

In order to prove convergence, we show a useful bound on $\dist_\W(\wb{w}_K):= \min_{w\in \W}\gnorm{w-\wb{w}_K}{}{}$.

Let $\wb{u}_K := \frac{1}{\Gamma_K} \tsum_{k = 0}^K\gamma_ku_{k+1}$. Then,

\begin{equation}\label{eq:int_rel_prox_23}
\dist_\W(\wb{w}_K)^2\le \gnorm{\wb{u}_K - \wb{w}_K}{}{2} = \frac{1}{\Gamma_K^2}\gnorm{\tsum_{k = 0}^K\gamma_k\boldsymbol{\xi}_{k+1}}{}{2}.
\end{equation}

Moreover, we have

\begin{equation}\label{eq:int_rel_prox_24}
\EE [\dist_\W(\wb{w}_K)^2] \le \frac{1}{\Gamma_K^2}\tsum_{k = 0}^K\gamma_k^2\EE\gnorm{\boldsymbol{\xi}_{k+1}}{}{2}
\end{equation}

Note that $\wb{u}_K \in \W$. Hence, first relation in (<ref>) follows by definition of $\dist_\W(\cdot)$ function. Equality follows from definition of $\wb{u}_K$ and $\wb{w}_K$.
To obtain (<ref>), note that $\{\boldsymbol{\xi}_k\}_{k = 1}^{K+1}$ are i.i.d. random variable with mean $0$. Expanding $\gnorm{\tsum_{k = 0}^K \gamma_k\boldsymbol{\xi}_{k+1}}{}{2}$, using linearity of expectation and noting $\EE[\boldsymbol{\xi}_i^T\boldsymbol{\xi}_j] = 0$ for all $i \ne j$, we conclude (<ref>). Hence, we conclude the proof.

We prove the following convergence rate result for Algorithm <ref> for the risk of SVI/SSP problem. In particular, we assume that the algorithm performs a single-pass over the dataset $\bS \sim \Pcal^n$ containing $n$ i.i.d. datapoints.

Consider the stochastic (<ref>) problem with operators $F_{\bbeta} \in \M^1(L, M)$. Let $\A$ be the single-pass NISPP method (Algorithm <ref>) where sequence $\{\gamma_k\}_{k \ge 0}, \{\lambda_k\}_{k \ge0}$ satisfy

\begin{equation}\label{eq:step_size_condition}
\gamma_k \lambda_k = \gamma_0\lambda_0
\end{equation}
for all $k \ge 0$. Moreover, $\bB_{k+1}$ are independent samples from a product distribution $\bB_{k+1} \sim \Pcal^{B_{k+1}}$ and $\bB_{k+1} \subset \bS$. Then, we have

\begin{align}
\VIgap{\A}{F} \le \nu + \frac{Z_0(K)}{\Gamma_K}  + M \sqrt{\frac{1}{\Gamma_K^2}\tsum_{k = 0}^K\gamma_k^2\sigma_{k+1}^2d}, \label{eq:conv_rate_svi_risk}
\end{align}

where, $Z_0(K) := \frac{3\gamma_0\lambda_0}{2}D^2 +  \frac{4M^2 + 3L^2D^2}{\gamma_0\lambda_0}\tsum_{k = 0}^K\gamma_k^2 + \frac{5\gamma_0\lambda_0d}{2}\tsum_{k =1}^K\sigma_{k+1}^2$ and $\Gamma_K := \tsum_{k = 0}^K \gamma_k$.

Similarly, $\A$ applied to stochastic (<ref>) problem achieves

\begin{equation}\label{eq:conv_spp_risk}
\SPgap{\A}{f} \le  \nu +  \frac{Z_0(K)}{\Gamma_K}  + M \sqrt{\frac{1}{\Gamma_K^2}\tsum_{k = 0}^K\gamma_k^2\sigma_{k+1}^2d}.
\end{equation}

Let $w\in\W$. Then

\begin{align}
\inprod{F(w)}{w_{k+1}-w} &= \inprod{F(w)}{u_{k+1}-w} + \inprod{F(w)}{w_{k+1}-u_{k+1}}.\label{eq:int_rel_prox_7}
\end{align}

We will analyze each term above separately. First, note that

\begin{align}
\inprod{F(w)}{w_{k+1}-u_{k+1}} &\le \frac{1}{2\lambda_k}\gnorm{F(w)}{}{2} + \frac{\lambda_k}{2}\gnorm{\boldsymbol{\xi}_{k+1}}{}{2} \le \frac{M^2}{2\lambda_k} + \frac{\lambda_k}{2}\gnorm{\boldsymbol{\xi}_{k+1}}{}{2}.\label{eq:int_rel_prox_6}
\end{align}

Note that

\begin{align}
\inprod{F(w)}{u_{k+1}-w} &\le \inprod{F(u_{k+1})}{u_{k+1}-w}\nonumber\\
&= \inprod{F_{\bB_{k+1}}(u_{k+1})}{u_{k+1}-w} + \inprod{F(u_{k+1}) - F_{\bB_{k+1}}(u_{k+1})}{u_{k+1}-w}\nonumber\\
&\le \lambda_k\inprod{u_{k+1}-w_k}{w - u_{k+1}} + \nu + \inprod{F(u_{k+1}) - F_{\bB_{k+1}}(u_{k+1})}{u_{k+1}-w}\nonumber\\
&=\frac{\lambda_k}{2} \bracket{\gnorm{w-w_k}{}{2} - \gnorm{w-u_{k+1}}{}{2} - \gnorm{u_{k+1} - w_k}{}{2}} + \nu \nonumber\\
&\qquad+ \inprod{F(u_{k+1}) - F_{\bB_{k+1}}(u_{k+1})}{u_{k+1}-w}\label{eq:int_rel_prox_2},
\end{align}

where first inequality follows from monotonicity and second inequality follows from (<ref>).
Now note that

\begin{align*}
&\inprod{F(u_{k+1}) - F_{\bB_{k+1}}(u_{k+1})}{u_{k+1}-w} \\
&\quad= \inprod{F(u_{k+1}) - F_{\bB_{k+1}}(u_{k+1})- [F(w_k) - F_{\bB_{k+1}}(w_{k})]}{u_{k+1}-w}\\
&\quad\quad + \inprod{F(w_k) - F_{\bB_{k+1}}(w_{k})}{u_{k+1} - w_k} +
\inprod{F(w_k) - F_{\bB_{k+1}}(w_{k})}{w_k-w}\\
&\quad \le \frac{\lambda_k}{6L^2}\gnorm{F(u_{k+1}) - F(w_k)}{}{2} + \frac{3L^2}{2\lambda_k}\gnorm{u_{k+1}-w}{}{2}\\
&\quad\quad + \frac{\lambda_k}{6L^2}\gnorm{F_{\bB_{k+1}}(u_{k+1}) - F_{\bB_{k+1}}(w_{k})}{}{2} + \frac{3L^2}{2\lambda_k}\gnorm{u_{k+1}-w}{}{2}\\
&\quad\quad+  \frac{3}{2\lambda_k}\gnorm{F(w_k) - F_{\bB_{k+1}}(w_{k})}{}{2} +  \frac{\lambda_k}{6} \gnorm{u_{k+1}- w_k}{}{2}+
\inprod{F(w_k) - F_{\bB_{k+1}}(w_{k})}{w_k-w}\\
&\quad \le \frac{\lambda_k}{2}\gnorm{u_{k+1}- w_k}{}{2} + \frac{3L^2}{\lambda_k}\gnorm{u_{k+1}-w}{}{2} + \frac{3}{2\lambda_k}\gnorm{F(w_k) - F_{\bB_{k+1}}(w_{k})}{}{2} \\
\inprod{F(w_k) - F_{\bB_{k+1}}(w_{k})}{w_k-w},
\end{align*}

where last inequality follows from $L$-Lipschitz continuity of $F$ and $F_{\bB_{k+1}}$.
Noting that $\gnorm{u_{k+1}-w}{}{} \le D$ for all $w \in \W$ and using the above bound in (<ref>), we have

\begin{align}
\inprod{F(w)}{u_{k+1}-w} &\le \inprod{F(u_{k+1})}{u_{k+1}-w} \nonumber\\
&\le \frac{\lambda_k}{2}\bracket{ \gnorm{w-w_k}{}{2} - \gnorm{w-u_{k+1}}{}{2} } \nonumber\\
&\quad+ \underbrace{\nu + \frac{3L^2D^2}{\lambda_k} + \frac{3}{2\lambda_k}\gnorm{F(w_k) - F_{\bB_{k+1}}(w_{k})}{}{2} +
\inprod{F(w_k) - F_{\bB_{k+1}}(w_{k})}{w_k-w}}_{E_k}
%&\le \frac{\lambda_k}{2}\bracket{ \gnorm{w-w_k}{}{2} - \gnorm{w-u_{k+1}}{}{2} } \underbrace{ + \nu + \frac{3(2L^2D^2 + M^2)}{2\lambda_k}+ \inprod{F(w_k) - F_{\bB_{k+1}}(w_{k})}{w_k-w} }_{E_k}
\label{eq:int_rel_prox_3},
\end{align}

Letting $u_0 := w_0$ and consequently $\boldsymbol{\xi}_0 = {\bf 0}$, we have from (<ref>)

\begin{align}
\inprod{F(w)}{u_{k+1}-w} &\le 	\inprod{F(u_{k+1})}{u_{k+1}-w}\nonumber\\  &\le \frac{\lambda_k}{2} \bracket{\gnorm{w-u_k}{}{2} - \gnorm{w-u_{k+1}}{}{2} + 2\inprod{w-u_k}{u_k-w_k} + \gnorm{u_k-w_k}{}{2}} + E_k \nonumber\\
&=\frac{\lambda_k}{2} \bracket{\gnorm{w-u_k}{}{2} - \gnorm{w-u_{k+1}}{}{2} } + \lambda_k \inprod{\boldsymbol{\xi}_k}{u_k-w} + \frac{1}{2}\lambda_k\gnorm{\boldsymbol{\xi}_k}{}{2} + E_k \label{eq:int_rel_prox_5}
\end{align}

Let us define an auxiliary sequence $\{z_k\}_{k \ge 0}$, where $z_0 = w_0$ and for all $ k \ge 1$, we have

\[z_k = \proj_{\W}[z_{k-1} - \boldsymbol{\xi}_k].\]

Then, due to the mirror-descent bound, we have

\begin{equation}
\tsum_{k =1}^K \inprod{\boldsymbol{\xi}_k}{z_k - w} \le \frac{1}{2}\gnorm{w-w_0}{}{2} - \frac{1}{2}\gnorm{w-z_K}{}{2} + \tsum_{k =1}^K \gnorm{\boldsymbol{\xi}_k}{}{2}.\label{eq:int_rel_prox_4}
\end{equation}

Moreover, noting that

\begin{align*}
\inprod{\boldsymbol{\xi}_k}{u_k - z_k}  &= \inprod{\boldsymbol{\xi}_k}{u_k - z_{k-1}}+\inprod{\boldsymbol{\xi}_k}{ z_{k-1} -z_k} \le \inprod{\boldsymbol{\xi}_k}{u_k - z_{k-1}} + \gnorm{\boldsymbol{\xi}_k}{}{}\gnorm{z_k - z_{k-1}}{}{}\\
&\le \inprod{\boldsymbol{\xi}_k}{u_k - z_{k-1}} + \gnorm{\boldsymbol{\xi}_k}{}{2}.
\end{align*}

Combining above relation with (<ref>), we have

\begin{equation}\label{eq:int_rel_prox_19}
\tsum_{k =1}^K\inprod{\boldsymbol{\xi}_k}{u_k - w} \le \frac{1}{2}\gnorm{w-w_0}{}{2} + 2\tsum_{k =1}^K\gnorm{\boldsymbol{\xi}_k}{}{2} + \tsum_{k =1}^K\inprod{\boldsymbol{\xi}_k}{u_k - z_{k-1}}.
\end{equation}

Now multiplying (<ref>) by $\gamma_k$ then summing from $k = 0$ to $K$; noting the definition of $\wb{w}_K$ and $\Gamma_K$; and using (<ref>), (<ref>) and (<ref>) along with assumption (<ref>) implies

\begin{align}
\Gamma_k\inprod{F(w)}{\wb{w}_K - w} &\le \gamma_0\lambda_0 \gnorm{w-w_0}{}{2} \nonumber\\
&\quad +\tsum_{k = 0}^K \bracket{\frac{\gamma_kM^2}{2\lambda_k} + \frac{\gamma_0\lambda_0}{2} \gnorm{\boldsymbol{\xi}_{k+1}}{}{2}  + \frac{5\gamma_0\lambda_0}{2} \gnorm{\boldsymbol{\xi}_k}{}{2} +\gamma_kE_k + \gamma_0\lambda_0\inprod{\boldsymbol{\xi}_k}{u_k - z_{k-1}}}\nonumber\\
\Rightarrow 	\Gamma_k\inprod{F(w)}{\proj_{\W}[\wb{w}_K] - w} &\le \gamma_0\lambda_0 \gnorm{w-w_0}{}{2} +\Gamma_K \inprod{F(w)}{\proj_{\W}[\wb{w}_K] - \wb{w}_K} \nonumber\\
&\quad +\tsum_{k = 0}^K \bracket{\frac{\gamma_kM^2}{2\lambda_k} + \frac{\gamma_0\lambda_0}{2} \gnorm{\boldsymbol{\xi}_{k+1}}{}{2}  + \frac{5\gamma_0\lambda_0}{2} \gnorm{\boldsymbol{\xi}_k}{}{2} +\gamma_kE_k + \gamma_0\lambda_0\inprod{\boldsymbol{\xi}_k}{u_k - z_{k-1}}}\nonumber\\
&\le \gamma_0\lambda_0 \gnorm{w-w_0}{}{2} + \Gamma_KM\dist_\W(\wb{w}_K) \nonumber\\
&\quad +\tsum_{k = 0}^K \bracket{\frac{\gamma_kM^2}{2\lambda_k} + \frac{\gamma_0\lambda_0}{2} \gnorm{\boldsymbol{\xi}_{k+1}}{}{2}  + \frac{5\gamma_0\lambda_0}{2} \gnorm{\boldsymbol{\xi}_k}{}{2} +\gamma_kE_k + \gamma_0\lambda_0\inprod{\boldsymbol{\xi}_k}{u_k - z_{k-1}}}\label{eq:int_rel_prox_12}
\end{align}

Now note that

\begin{align}
\tsum_{k = 0}^K \gamma_kE_k &= \Gamma_K \nu + \frac{3L^2D^2}{\gamma_0\lambda_0}\tsum_{k = 0}^K \gamma_k^2 + \frac{3}{2\gamma_0\lambda_0}\tsum_{k = 0}^K \gamma_k^2\gnorm{F_{\bB_{k+1}}(w_{k})-F(w_k)}{}{2}\nonumber\\
&\quad+ \gamma_0\lambda_0 \tsum_{k = 0}^K \inprod{\frac{1}{\lambda_k}(F(w_k) - F_{\bB_{k+1}}(w_{k}))}{w_k-w}\label{eq:int_rel_prox_9}
\end{align}

Define $\boldsymbol{\Delta}_k := \frac{1}{\lambda_k}(F(w_k) - F_{\bB_{k+1}}(w_{k}))$. Note that $\EE_{\bB_{k+1}}[\boldsymbol{\Delta}_k| w_k] = 0$. Moreover, define an auxiliary sequence $\{h_k\}_{k \ge 0}$ with $h_0 := w_0$ and

\[h_{k+1} := \proj_{\W}[h_k - \boldsymbol{\Delta}_k].\]
Then due to mirror descent bound, we have

\begin{align}
\tsum_{k = 0}^K \inprod{\boldsymbol{\Delta}_k}{h_{k+1}-w} \le \frac{1}{2}\gnorm{w-w_0}{}{2} - \frac{1}{2}\gnorm{w-h_{K+1}}{}{2} + \tsum_{k = 0}^K\gnorm{\boldsymbol{\Delta}_k}{}{2}.\label{eq:int_rel_prox_8}
\end{align}

\begin{align*}
\inprod{\boldsymbol{\Delta}_k}{h_k-h_{k+1}} \le \gnorm{\boldsymbol{\Delta}_k}{}{}\gnorm{h_k-h_{k+1}}{}{} \le \gnorm{\boldsymbol{\Delta}_k}{}{2}.
\end{align*}

Using above relation along with (<ref>), we have

\begin{align*}
\tsum_{k = 0}^K\inprod{\boldsymbol{\Delta}_k}{h_k - w} &= \tsum_{k = 0}^K[ \inprod{\boldsymbol{\Delta}_k}{h_{k+1} - w} + \inprod{\boldsymbol{\Delta}_k}{h_k - h_{k+1}}]\\
&\le \frac{1}{2}\gnorm{w-w_0}{}{2} - \frac{1}{2}\gnorm{w-h_{K+1}}{}{2} + 2\tsum_{k = 0}^K\gnorm{\boldsymbol{\Delta}_k}{}{2}\\
\Rightarrow\tsum_{k = 0}^K\inprod{\boldsymbol{\Delta}_k}{w_k - w}&= \tsum_{k = 0}^K[ \inprod{\boldsymbol{\Delta}_k}{w_k-h_k} + \inprod{\boldsymbol{\Delta}_k}{h_k-w}]\\
&\le \tsum_{k=0}^K\inprod{\boldsymbol{\Delta}_k}{w_k-h_k} + \frac{1}{2}\gnorm{w-w_0}{}{2} - \frac{1}{2}\gnorm{w-h_{K+1}}{}{2} + 2\tsum_{k = 0}^K\gnorm{\boldsymbol{\Delta}_k}{}{2}.
\end{align*}

Using the above relation in (<ref>), we have

\begin{align}
\tsum_{k = 0}^K \gamma_kE_k &= \Gamma_K \nu + \frac{3L^2D^2}{\gamma_0\lambda_0}\tsum_{k = 0}^K \gamma_k^2 + \frac{3}{2\gamma_0\lambda_0}\tsum_{k = 0}^K \gamma_k^2\gnorm{F_{\bB_{k+1}}(w_{k})-F(w_k)}{}{2}\nonumber\\
&\quad+\gamma_0\lambda_0 \tsum_{k =0}^K \inprod{\boldsymbol{\Delta}_k}{w_k-w} \nonumber\\
&\le \frac{\gamma_0\lambda_0}{2}\gnorm{w-w_0}{}{2} + \Gamma_K\nu \frac{3L^2D^2}{\gamma_0\lambda_0}\tsum_{k = 0}^K \gamma_k^2 + \frac{3}{2\gamma_0\lambda_0}\tsum_{k = 0}^K \gamma_k^2\gnorm{F_{\bB_{k+1}}(w_{k})-F(w_k)}{}{2}\nonumber\\
&\quad+ 2\gamma_0\lambda_0\tsum_{k = 0}^{K} \gnorm{\boldsymbol{\Delta}_k}{}{2} + \gamma_0\lambda_0\tsum_{k = 0}^K \inprod{\boldsymbol{\Delta}_k}{w_k-h_k}.\label{eq:int_rel_prox_13}
\end{align}

Finally, note that for all valid $k$, we have

\begin{align}
\EE[\gnorm{\bg_k}{}{2}] &= \sigma_k^2d\label{eq:int_rel_prox_10},\\
\EE[\gnorm{\boldsymbol{\Delta}_k}{}{2}] &=\EE_{w_k}[\EE_{\bB_{k+1}}[\gnorm{\boldsymbol{\Delta}_k}{}{2}| w_k]] \nonumber\\
&= \frac{1}{\lambda_k^2}\EE_{w_k} \EE_{\bB_{k+1}}\bracket{\gnorm{F_{\bB_{k+1}}(w_{k}) - F(w_k)}{}{2} | w_k}\nonumber\\
&\le \frac{1}{\lambda_k^2}\EE_{w_k}\EE_{\bB_{k+1}}\gnorm{F_{\bB_{k+1}}(w_{k})}{}{2} \le \frac{M^2}{\lambda_k^2}\label{eq:int_rel_prox_11},\\
\EE[\inprod{\boldsymbol{\Delta}_k}{w_k-h_k}] &= \EE[\inprod{\EE[\boldsymbol{\Delta}_k|w_k, h_k]}{w_k-h_k}] = 0 \label{eq:int_rel_prox_14},\\
\EE[\inprod{\boldsymbol{\xi}_k}{u_k-z_{k-1}}] &= \EE[\inprod{\EE[\boldsymbol{\xi}_k| u_k, z_{k-1}]}{u_k-z_{k-1}}] = 0\label{eq:int_rel_prox_15},
\end{align}

where, in (<ref>), we used the fact that $F_{\bbeta})$ is $M$-bounded for all $\boldsymbol{\beta} \in \bS$.

Now, using (<ref>) in relation (<ref>), noting the bound on $\dist_\W(\wb{w}_K)$ from Proposition <ref> (in particular (<ref>)), taking supremum with respect to $w \in \W$, then taking expectation and noting (<ref>)-(<ref>), we have

\begin{align*}
\Gamma_k\EE\sup_{w \in \W}\inprod{F(w)}{\proj_{\W}[\wb{w}_K] - w} &\le \frac{3\gamma_0\lambda_0}{2} D^2 + \tsum_{k =0}^K\bracket{\frac{\gamma_k(4M^2+3L^2D^2)}{\lambda_k} + \frac{5\gamma_0\lambda_0}{2}\sigma_{k+1}^2d+ \gamma_k\nu} \\
&\quad+ \Gamma_K M \sqrt{\frac{1}{\Gamma_K^2}\tsum_{k = 0}^K\gamma_k^2\sigma_{k+1}^2d}.\\
\Rightarrow \EE\sup_{w \in \W}\inprod{F(w)}{\proj_{\W}[\wb{w}_K] - w} &\le \nu +  \frac{1}{\Gamma_K} [\frac{3\gamma_0\lambda_0}{2}D^2 +  \frac{4M^2 + 3L^2D^2}{\gamma_0\lambda_0}\tsum_{k = 0}^K\gamma_k^2 + \frac{5\gamma_0\lambda_0d}{2}\tsum_{k =1}^K\sigma_{k+1}^2 ]\\
&\quad + M \sqrt{\frac{1}{\Gamma_K^2}\tsum_{k = 0}^K\gamma_k^2\sigma_{k+1}^2d},
\end{align*}

where in the first inequality, we used the fact that $\EE[\dist_\W(\wb{w}_K)] \le \sqrt{\EE[\dist_\W(\wb{w}_K)^2]}$. Hence, we conclude the proof of (<ref>).

Now, we extend this for (<ref>). Denote $u^{k+1} = (\wt{x}^{k+1}, \wt{y}^{k+1})$. Then, we have

\begin{align*}
\inprod{F(u_{k+1})}{u_{k+1} -w}&= \inprod{\grad_x f(\wt{x}_{k+1}, \wt{y}_{k+1} )}{\wt{x}_{k+1} - x} + \inprod{-\grad_y f(\wt{x}_{k+1}, \wt{y}_{k+1} )}{\wt{y}_{k+1} - y}\\
&\ge f(\wt{x}_{k+1}, \wt{y}_{k+1}) - f(x, \wt{y}_{k+1})  + [-f(\wt{x}_{k+1}, \wt{y}_{k+1}) + f(\wt{x}_{k+1}, y)]\\
&= f(\wt{x}_{k+1}, y) - f(x, \wt{y}_{k+1}).
\end{align*}

Using the above in (<ref>), we obtain,

\begin{align*}
f(\wt{x}_{k+1}, y) - f(x, \wt{y}_{k+1}) &\le \frac{\lambda_k}{2} \bracket{\gnorm{w-u_k}{}{2} - \gnorm{w-u_{k+1}}{}{2} } + \lambda_k \inprod{\boldsymbol{\xi}_k}{u_k-w} + \frac{1}{2}\lambda_k\gnorm{\boldsymbol{\xi}_k}{}{2} + E_k.
%	\\
%	\Rightarrow
%	f(x_{k+1}, y; \bS) -f(x, y_{k+1}; \bS) &\le \frac{\lambda_k}{2} \bracket{\gnorm{w-\wt{w}_k}{}{2} - \gnorm{w-u_{k+1}}{}{2} } + \lambda_k \inprod{\boldsymbol{\xi}_k}{\wt{w}_k-w} + \frac{1}{2}\lambda_k\gnorm{\boldsymbol{\xi}_k}{}{2} + E_k\\
%	&\qquad + [f(x_{k+1}, y; \bS) - f(\wt{x}_{k+1}, y; \bS)] - [f(x, y_{k+1}; \bS) - f(x, \wt{y}_{k+1}; \bS)]\\
%	&\le \frac{\lambda_k}{2} \bracket{\gnorm{w-\wt{w}_k}{}{2} - \gnorm{w-u_{k+1}}{}{2} } + \lambda_k \inprod{\boldsymbol{\xi}_k}{\wt{w}_k-w} + \frac{1}{2}\lambda_k\gnorm{\boldsymbol{\xi}_k}{}{2} + E_k\\
\end{align*}

Now, using Proposition <ref> to bound the distance between points $\tfrac{1}{\Gamma_K}(\tsum_{k =0}^K\gamma_k\wt{x}_{k+1}, \tsum_{k =0}^K\gamma_k\wt{y}_{k+1})$ and $(\Pi_\X[\wb{x}_K], \Pi_\Y[\wb{y}_K])$ and using Jensen's inequality to conclude that

\begin{align*}
\frac{1}{\Gamma_K}\tsum_{k = 0}^K\gamma_k [f(\wt{x}_{k+1}, y) - f(x, \wt{y}_{k+1})] &\ge f(\tfrac{1}{\Gamma_K}\tsum_{k = 0}^K\gamma_k \wt{x}_{k+1}, y) - f(x, \tfrac{1}{\Gamma_K} \tsum_{k = 0}^K\gamma_k\wt{y}_{k+1})]\\
&\ge f(\wb{x}_K,y) - f(x, \wb{y}_K) - M\gnorm*{\begin{bmatrix}
\tfrac{1}{\Gamma_K}\tsum_{k = 0}^K\gamma_k \wt{x}_{k+1} - \proj_\X[\wb{x}_K]\\
\tfrac{1}{\Gamma_K}\tsum_{k = 0}^K\gamma_k \wt{y}_{k+1} - \proj_\Y[\wb{y}_K]
\end{bmatrix}}{}{}
\end{align*}
and retracing the steps of this proof from (<ref>), we obtain (<ref>). Hence, we conclude the proof.

§.§ Differential privacy of the NISPP method

First, we show a simple bound on $\ell_2$-sensitivity for updates of NISPP method.

Suppose $\nu \le \frac{2M^2}{\lambda_kB_{k+1}^2}$ then $\ell_2$-sensitivity of updates of Algorithm <ref>
is at most $\frac{4M}{\lambda_kB_{k+1}}$

where $B_{k+1} = \abs{\bB_{k+1}}$ is the batch size of $k$-th iteration.

Let $w_k$ be an iterate in the start of $k$-th iteration of Algorithm <ref>. Suppose $\bB_{k+1}$ and $\bB_{k+1}'$ be two different batches used in $k$-th iteration to obtain $u_{k+1}$ and $u_{k+1}'$, respectively. Also note that $\bB_{k+1}$ and $\bB_{k+1}'$ differ in only single datapoint. Then, due to (<ref>), we have for all $w \in \W$

\begin{align*}
\inprod{F_{\bB_{k+1}}(u_{k+1}) + \lambda_k(u_{k+1} - w_k) }{w - u_{k+1}} &\ge -\nu\\
\inprod{F_{\bB_{k+1}'}(u_{k+1}') + \lambda_k(u_{k+1}' - w_k) }{w - u_{k+1}'} &\ge -\nu
\end{align*}

Using $w = u_{k+1}'$ in the first relation and $w = u_{k+1}$ in the second relation above and then summing, we obtain

\begin{align*}
\inprod{ F_{\bB_{k+1}}(u_{k+1}) -F_{\bB_{k+1}'}(u_{k+1}') }{ u_{k+1} - u_{k+1}'} &\le 2\nu - \lambda_k \gnorm{u_{k+1} - u_{k+1}'}{}{2}\\
\Rightarrow \inprod{F_{\bB_{k+1}}(u_{k+1}) - F_{\bB_{k+1}}(u_{k+1}') }{ u_{k+1} - u_{k+1}'}&\le \inprod{  F_{\bB_{k+1}'}(u_{k+1}')- F_{\bB_{k+1}}(u_{k+1}') }{ u_{k+1} - u_{k+1}'} \\
&\quad+ 2\nu - \lambda_k \gnorm{u_{k+1} - u_{k+1}'}{}{2}
\end{align*}

Now, noting that $F_{\bB_{k+1}}$ is a monotone operator and denoting $\mathbf{a}_{k+1} := \gnorm{ F_{\bB_{k+1}'}(u_{k+1}')- F_{\bB_{k+1}}(u_{k+1}')}{}{}$, $p_{k+1} := \gnorm{w_{k+1} - w'_{k+1}}{}{} = \gnorm{u_{k+1} - u'_{k+1}}{}{}$ we have

\begin{align}
0 &\le \inprod{  F_{\bB_{k+1}'}(u_{k+1}')- F_{\bB_{k+1}}(u_{k+1}') }{ u_{k+1} - u_{k+1}'}+ 2\nu - \lambda_k \gnorm{u_{k+1} - u_{k+1}'}{}{2}\nonumber\\
&\le \mathbf{a}_{k+1} p_{k+1} - \lambda_k p_{k+1}^2 + 2\nu. \label{eq:int_rel_prox_16}
\end{align}

Finally noting that if $\boldsymbol{\beta}$ and $\boldsymbol{\beta}'$ are the differing datapoints in $\bB_{k+1}$ and $\bB_{k+1}'$, then

\begin{align*}
\mathbf{a}_{k+1} &= \frac{1}{B_{k+1}}\gnorm{F_{\bbeta'}(u'_{k+1}) - F_{\bbeta}(u'_{k+1}) }{}{} \le \frac{2M}{B_{k+1}}.
\end{align*}

Using the above relation in (<ref>) and noting that $\ell_2$-sensitivity $p_{k+1} = \gnorm{w_{k+1}-w_{k+1}'}{}{} = \gnorm{u_{k+1} -u_{k+1}'}{}{}$, we have, $p_{k+1}$ satisfies

\[p_{k+1}^2 -\frac{2M}{\lambda_k B_{k+1}} p_{k+1} - \frac{2\nu}{\lambda_k} \le 0.\]

This implies

\[p_{k+1} \le \frac{M}{\lambda_kB_{k+1}} + \sqrt{\frac{M^2}{\lambda_k^2B_{k+1}^2} + \frac{2\nu}{\lambda_k}} \le \frac{2M}{\lambda_kB_{k+1}} + \sqrt{\frac{2\nu}{\lambda_k}}.\]

Setting $\nu \le \frac{2M^2}{\lambda_kB_{k+1}^2}$, we have $p_{k+1} \le \frac{4M}{\lambda_kB_{k+1}}.$
Hence, we conclude the proof.

Using the $\ell_2$-sensitivity result above along with Proposition <ref> and <ref>, we immediately obtain the following:

Algorithm <ref> with batch sizes $(B_{k+1})_{k \in [K]_0}$, parameters $(\lambda_k)_{k \in [K]_0}$, variance $\sigma_{k+1}^2 = \frac{32M^2}{\lambda_k^2B_{k+1}^2}\frac{\ln(1/\eta)}{\eps^2}$ and $\nu$ satisfying assumptions of Proposition <ref> is $(\eps, \eta)$-differentially private.

Note that for any $k \ge 0$, we have

\begin{align}
\exp\braces*{(\alpha-1) D_\alpha(w_{k+1} || w_{k+1}')} &= \EE_{\mu_{k+1}'}\bracket*{\paran*{\frac{d\mu_{k+1}}{d\mu_{k+1}'}}^{\alpha}} \nonumber \\
&= \tsum_{t = 0}^{K} \PP(i \in \bB_{t+1})\EE_{\mu_{k+1}'(\cdot| i \in \bB_{t+1})} \bracket*{\paran*{\frac{d\mu_{k+1}(\cdot| i \in \bB_{t+1})}{d \mu_{k+1}'( \cdot| i \in \bB_{t+1})}}^\alpha } \nonumber \\
&= \tsum_{t = 0}^{k} \PP(i \in \bB_{t+1})\EE_{\mu_{k+1}'(\cdot| i \in \bB_{t+1})} \bracket*{\paran*{\frac{d\mu_{k+1}(\cdot| i \in \bB_{t+1})}{d \mu_{k+1}'( \cdot| i \in \bB_{t+1})}}^\alpha } \nonumber \\
&\quad +  \tsum_{t = k+1}^{K} \PP(i \in \bB_{t+1})\EE_{\mu_{k+1}'(\cdot| i \in \bB_{t+1})} \bracket*{\paran*{\frac{d\mu_{k+1}(\cdot| i \in \bB_{t+1})}{d \mu_{k+1}'( \cdot| i \in \bB_{t+1})}}^\alpha }, \label{eq:int_rel_prox_17}
\end{align}

where second equality follows from law of total expectation.
Now note that conditioned on $i \in \bB_{t+1}$ for $t \ge k+1$ then $\mu_{k+1}$ and $\mu_{k+1}'$ have the same distribution. Hence,  we have

\begin{align}\label{eq:int_rel_prox_18}
&\EE_{\mu_{k+1}'(\cdot| i \in \bB_{t+1})} \bracket*{\paran*{\frac{d\mu_{k+1}(\cdot| i \in \bB_{t+1})}{d \mu_{k+1}'( \cdot| i \in \bB_{t+1})}}^\alpha }= 1 & \forall \ t \ge k+1.
\end{align}

On the other hand, by inductively applying the data-processing inequality for the Rényi divergence, we have

\begin{align}\label{eq:int_rel_prox_20}
&\EE_{\mu_{k+1}'(\cdot| i \in \bB_{t+1})} \bracket*{\paran*{\frac{d\mu_{k+1}(\cdot| i \in \bB_{t+1})}{d \mu_{k+1}'( \cdot| i \in \bB_{t+1})}}^\alpha } \le \EE_{\mu_{t+1}'(\cdot| i \in \bB_{t+1})} \bracket*{\paran*{\frac{d\mu_{t+1}(\cdot| i \in \bB_{t+1})}{d \mu_{t+1}'( \cdot| i \in \bB_{t+1})}}^\alpha } &\forall \ t \le k.
\end{align}

Next, let $C_0 = 0$, $C_{k+1} = \sum_{t = 0}^k B_{t+1}$ and note that $\bB_{t+1}$ is sampled without replacement (i.e., case of single-pass algorithm), we have

\begin{align}\label{eq:int_rel_prox_21}
\PP(i \in \bB_{t+1}) = \frac{B_{t+1}}{n -C_t} \prod_{s = 1}^{t} \frac{n- C_s}{n- C_{s-1}} = \frac{B_{t+1}}{n} .
\end{align}

Using (<ref>), (<ref>) and (<ref>) into (<ref>), we have

\begin{align*}
\exp \braces*{(\alpha-1)D_\alpha(w_{k+1}|| w_{k+1}')} \le \tsum_{t= 0}^k \frac{B_{t+1}}{n} \EE_{\mu_{t+1}'(\cdot| i \in \bB_{t+1})} \bracket*{\paran*{\frac{d\mu_{t+1}(\cdot| i \in \bB_{t+1})}{d \mu_{t+1}'( \cdot| i \in \bB_{t+1})}}^\alpha } + \tsum_{t = k+1}^K\frac{B_{t+1}}{n}
\end{align*}

Now note that from the Rényi divergence bound for Gaussians with different mean MISSING REF,

$$ \EE_{\mu_{t+1}'(\cdot| i \in \bB_{t+1})} \bracket*{\paran*{\frac{d\mu_{t+1}(\cdot| i \in \bB_{t+1})}{d \mu_{t+1}'( \cdot| i \in \bB_{t+1})}}^\alpha }\leq \exp\Big\{(\alpha-1) \frac{\alpha p_{t+1}^2}{2\sigma_{t+1}^2}\Big\},$$

where $p_{t+1}$ is the $\ell_2$-sensitivity of $w_{t+1}$. Hence, by Proposition <ref>

\begin{eqnarray*}
\EE_{\mu_{t+1}'(\cdot| i \in \bB_{t+1})} \bracket*{\paran*{\frac{d\mu_{t+1}(\cdot| i \in \bB_{t+1})}{d \mu_{t+1}'( \cdot| i \in \bB_{t+1})}}^\alpha }
&\leq& \exp\braces*{(\alpha-1) \frac{8\alpha M^2}{\lambda_t^2B_{t+1}^2\sigma_{t+1}^2}} = \exp\braces*{(\alpha-1) \frac{\alpha \rho^2}{2}},
\end{eqnarray*}

where in the last inequality we used $\sigma_{t+1} = \frac{4M}{\lambda_t B_{t+1}\rho}$.
Therefore, for any $\alpha > 1$, we conclude that

\begin{align*}
\exp \braces*{(\alpha-1) D_\alpha(w_{k+1}|| w_{k+1}')} &\le \tsum_{t = 0}^k \frac{B_{t+1}}{n}\exp\braces*{(\alpha-1)\frac{\alpha\rho^2}{2}} + \tsum_{t = k+1}^K \frac{B_{t+1}}{n} \\
&\le \exp\braces*{(\alpha-1)\frac{\alpha\rho^2}{2}} \paran*{\tsum_{t = 0}^K \frac{B_{t+1}}{n} }= \exp\braces*{(\alpha-1)\frac{\alpha\rho^2}{2}}
\end{align*}

where in the last relation, we used the fact that this is a single pass algorithm hence $\tsum_{t= 0}^KB_{t+1} = n$. Hence, we obtain for arbitrary $k \in \{0, \dots, K\}$ we have
$D_\alpha(w_{k+1}|| w_{k+1}') \le \frac{1}{2}\alpha\rho^2$. Using quasi-convexity of the Rényi divergence,
we have What is $[K]_0$ here?
$D_\alpha(\wb{w}_K || \wb{w}_K') \le \max_{k\in [K]_0} D_\alpha(w_{k+1}||w_{k+1}') \le \frac{1}{2}\alpha\rho^2$. Finally, we
we use again the data-processing inequality to conclude that $D_\alpha(\A(\bS) || \A(\bS')) \le \frac{1}{2}\alpha\rho^2$, which proves
the result.

Now, we provide a policy for setting $\gamma_k,  \lambda_k$ and $B_{k+1}$ to obtain population risk bounds for DP-SVI and DP-SSP problem by the NISPP method.

Algorithm <ref> with disjoint batches $\bB_{k+1}$ is of size $B_{k+1} = B := n^{1/3}$ for all $k \ge 0$ and the following parameters

γ_k    = 1,
λ_k = λ_0 := max{M/D, L}maxn^1/3, √(dln(1/η))/,

σ_k+1^2    = 32M^2/ Bλ_0 ln(1/η)/^2,
ν= 2M^2/λ_0B^2,

is  $(\eps, \eta)$-differentially private and achieves expected SVI-gap (SSP-gap, respectively)

\[%\EE_{\A, \bS} \big[\VIgap{\A(\bS)}{F}\big] \le
O\paran*{(M+LD)D\bracket*{\frac{1}{n^{1/3}} + \frac{\sqrt{d\ln(1/\eta)}}{\eps n^{2/3}}} %+\frac{LD^2}{n^{1/3}}
\]

Note that values of $\nu$, $\sigma_{k+1}$ and other required conditions proposed in Propositions <ref> and <ref> are satisfied. Hence, this algorithm is $(\eps, \eta)$-differentially private.

Moreover, all requirements of Theorem <ref> are satisfied. In order to maintain single pass over the dataset, we require $K = \frac{n}{B} = n^{2/3}$ iterations. Then,  we provide individual bounds on the terms of (<ref>) ((<ref>), respectively) and conclude the corollary using Theorem <ref>.

Note that we are using a constant parameter policy. Hence, $\sigma_{k+1} = \sigma = \frac{4M}{\rho B\lambda_0}$ for all $k \ge 0$. Substituting appropriate parameter values, we have

\begin{align*}
\nu &= \frac{2MD}{n^{2/3}\max\{n^{1/3}, \sqrt{d\ln(1/\eta)}/\eps\} } \le \frac{2MD}{n},\\
M \sqrt{\frac{1}{\Gamma_K^2}\tsum_{k = 0}^K\gamma_k^2\sigma_{k+1}^2d}&= \frac{M\sqrt{d}\sigma_{k+1}}{\sqrt{K}} = \frac{4M^2\sqrt{2d\ln(1/\eta)}}{\eps n^{2/3}\lambda_0} \le \frac{4\sqrt{2}MD}{n^{2/3}}, \\
\frac{3\lambda_0D^2}{2K} &\le \frac{3(M+LD)D}{2} \paran*{\frac{1}{n^{1/3}} + \frac{\sqrt{d\ln(1/\eta)}}{\eps n^{2/3}}},\\
\frac{4M^2 + 3L^2D^2}{\lambda_0} &\le %\frac{7M^2}{\lambda_0} \le
\frac{4MD}{n^{1/3}} + \frac{3LD^2}{n^{1/3}},\\
\frac{5\lambda_0d\sigma^2}{2} &= \frac{40M^2d\ln(1/\eta)}{\eps^2B^2 \lambda_0} \le \frac{40MD \sqrt{d\ln(1/\eta)}}{\eps n^{2/3}}.
\end{align*}

Substituting these bounds in Theorem <ref>, we conclude the proof.

We have the following remarks for NISPP method:

* In order to obtain $\nu$-approximate solution of the subproblem of NISPP method satisfying (<ref>), we can use the Operator Extrapolation (OE) method (see Theorem 2.3 [25]). OE method outputs a solution $u_{k+1}$ satisfying $\gnorm{u_{k+1} -w^*_{k+1}}{}{} \le \zeta$ in $ \frac{L+\lambda_0}{\lambda_0}\ln(\frac{D}{\zeta})$ iterations, where $w_{k+1}^*$ is an exact SVI solution for problem (<ref>). Furthermore, we have for all $w \in \W$,

\begin{align*}
0 &\le \inprod{F(w^*_{k+1}) + \lambda_k(w^*_{k+1}-w_k)}{w-w^*_{k+1}} \\
&= \inprod{F(u_{k+1}) + F(w^*_{k+1}) - F(u_{k+1}) + \lambda_k (u_{k+1}- w_k) + \lambda_k(w^*_{k+1}- u_{k+1}) }{w -w^*_{k+1}} \\
&\le \inprod{F(u_{k+1}) + \lambda_k (u_{k+1}- w_k)}{w -w^*_{k+1}} + (L+\lambda_k)\gnorm{u_{k+1}-w^*_{k+1}}{}{}\gnorm{w -w^*_{k+1}}{}{}\\
&\le  \inprod{F(u_{k+1}) + \lambda_k (u_{k+1}- w_k)}{w -w^*_{k+1}} + (L+\lambda_k)D\gnorm{u_{k+1}-w^*_{k+1}}{}{}\\
&=  \inprod{F(u_{k+1}) + \lambda_k (u_{k+1}- w_k)}{w- u_{k+1} + u_{k+1}-w^*_{k+1}} + (L+\lambda_k)D\gnorm{u_{k+1}-w^*_{k+1}}{}{}\\
&\le \inprod{F(u_{k+1}) + \lambda_k (u_{k+1}- w_k)}{w- u_{k+1}} + \gnorm{F(u_{k+1}) + \lambda_k (u_{k+1}- w_k)}{}{}\gnorm{u_{k+1}-w^*_{k+1}}{}{}\\
&\quad + (L+\lambda_k)D\gnorm{u_{k+1}-w^*_{k+1}}{}{}\\
&\le \inprod{F(u_{k+1}) + \lambda_k (u_{k+1}- w_k)}{w- u_{k+1}} + [LD+M+2\lambda_kD]\gnorm{u_{k+1}-w^*_{k+1}}{}{}
\end{align*}

$\zeta=\nu/[LD+M+2\lambda_kD]$, we obtain that $u_{k+1}$ is a $\nu$-approximate solution satisfying (<ref>). Using the convergence rate above, we require $\frac{L+\lambda_0}{\lambda_0}\ln\frac{MD+LD^2+2\lambda_kD^2}{\nu}$ operator evaluations.

Note that since, $\lambda_0 \ge L$,
we have $\frac{L+\lambda_0}{\lambda_0} \le 2$.

\begin{align}
\ln{\frac{MD+LD^2 + 2\lambda_kD^2}{\nu}} &\le \ln{\frac{4\lambda_kD^2}{\nu}}\label{eq:int_rel_prox_29}\\
&= \ln\paran[\big]{\frac{2\lambda_0^2D^2B^2}{M^2}} \nonumber\\
&= \ln\paran[\big]{ n^{2/3}\max\braces[\big]{n^{2/3}, \frac{d\ln(1\privb)}{\priva^2}}\max\braces[\big]{1, \frac{L^2D^2}{M^2}}  }\nonumber
\end{align}

Hence, each iteration of NISPP method requires $O(\log{n})$ iterations of OE method for solving the subproblem. Moreover, each iteration of the OE method requires $2B$ stochastic operator evaluations. Hence, we require $O(KB\log{n})$ stochastic operator evaluations in the entire run of NISPP (Algorithm <ref>). Noting that $KB = n$, we conclude that this is a near linear time algorithm and also performs only a single pass over the data in the stochastic outer-loop. We provide the details of OE method in the Appendix <ref>.

* For non-DP version of NISPP method, i.e., $\sigma_k = 0$ for all $k$, we can easily obtain population risk bound of $O(\frac{MD}{\sqrt{n}})$ by setting $\lambda_0 = \frac{M}{D}\sqrt{n}$, $B = 1$ (or $K = n$) and $\nu = \frac{MD}{\sqrt{n}}$ in Corollary <ref>.

In view of Corollary <ref>, it seems that running NISPP method for $n^{3/2}$ stochastic operator evaluations may provide optimal risk bounds. However, running that many stochastic operator evaluations requires multi-pass over the dataset so, in principle, this would only provide bounds in the empirical risk.

In order to compute the population risk of this multi-pass version, we analyze the stability of NISPP and provide generalization guarantees which result in optimal population risk.

§ STABILITY OF NISPP AND OPTIMAL RISK FOR DP-SVI  AND DP-SSP

In this section, we develop a multi-pass variant of NISPP method, and prove its stability to extrapolate empirical performance to population risk bounds.

§.§ Stability of NISPP method

Let us start with two adjacent datasets $\bS \simeq \bS'$. Suppose we run NISPP method on both datasets starting from the same point $w_0 \in \W$.
Then, in the following lemma, we provide bound on the how far apart trajectories of these two runs can drift.