Differentially Private Deep Learning with ModelMix

Training large neural networks with meaningful/usable differential privacy security guarantees is a demanding challenge.
In this paper, we tackle this problem by revisiting the two key operations in Differentially Private Stochastic Gradient Descent (DP-SGD): 1) iterative perturbation and 2) gradient clipping.
We propose a generic optimization framework, called ModelMix, which performs random aggregation of intermediate model states.
It strengthens the composite privacy analysis utilizing the entropy of the training trajectory
and improves the $(\epsilon, \delta)$ DP security parameters by an order of magnitude.

We provide rigorous analyses for both the utility guarantees and privacy amplification of ModelMix.
In particular, we present a formal study on the effect of gradient clipping in DP-SGD, which provides theoretical instruction on how hyper-parameters should be selected.
We also introduce a refined gradient clipping method, which can further sharpen the privacy loss in private learning when combined with ModelMix.

Thorough experiments with significant privacy/utility improvement are presented to support our theory.
We train a Resnet-20 network on CIFAR10 with $70.4\%$ accuracy via ModelMix given $(\epsilon=8, \delta=10^{-5})$ DP-budget, compared to the same performance but with $(\epsilon=145.8,\delta=10^{-5})$ using regular DP-SGD; assisted with additional public low-dimensional gradient embedding, one can further improve the accuracy to $79.1\%$ with $(\epsilon=6.1, \delta=10^{-5})$ DP-budget, compared to the same performance but with $(\epsilon=111.2, \delta=10^{-5})$ without ModelMix.

Differential Privacy; Rényi Differential Privacy;  Clipped Stochastic Gradient Descent; Deep Learning;

§ INTRODUCTION

Privacy concerns when learning with sensitive data are receiving increasing attention.
Many practical attacks have shown that without proper protection, the  model's parameters <cit.>, leakage on gradients during training <cit.>, or just observations on the prediction results <cit.> may enable an adversary to successfully distinguish and even reconstruct the private samples used for learning.
As an emergent canonical definition, Differential Privacy (DP) <cit.> provides a semantic privacy metric to quantify how hard it is for an adversary to infer the participation of an individual in an aggregate statistic.
As one of the most popular approaches, Differentially-Private Stochastic Gradient Descent (DP-SGD) <cit.> and its variants <cit.> have been widely studied over the last decade.
DP-SGD can be applied to almost all optimization problems in machine learning to produce rigorous DP guarantees without additional assumptions regarding the objective function or dataset.
However, despite its broad applicability, DP-SGD also suffers notoriously large utility loss especially when training cutting-edge deep models.
Its practical implementation is also known to be sensitive to hyper-parameter selections <cit.>. Indeed, even in theory, the effects of the two artificial privatization operations applied in DP-SGD, iterative gradient perturbation and gradient clipping, are still not fully-understood.

To understand why these two artificial modifications are the key to differentially privatize iterative methods, we need to first introduce the concept of sensitivity, which plays a key role in DP.
The sensitivity captures the maximum impact an individual sample from an input dataset may have on an algorithm's output.
It is the foundation of almost all DP mechanisms, including the Laplace/Gaussian and Exponential Mechanisms <cit.>, where the sensitivity determines how much randomization is needed to to hide any individual amongst the population with desired privacy.
Unfortunately, in many practical optimization problems, the end-to-end sensitivity is intractable or can only be loosely bounded.

To this end, DP-SGD proposes an alternative solution by assuming a more powerful adversary.
In most private (centralized) learning applications, the standard black-box adversary can only observe the final model revealed.
DP-SGD, on the other hand, assumes an adversary who can observe the intermediate updates during training. For convenience, we will call such an adversary a white-box adversary.

Provided such an empowered adversary, DP-SGD clips the gradient evaluated by each individual sample and adds random noises to the updates in each iteration.
Clipping guarantees that the sensitivity is bounded within each iteration.
The total privacy loss is then upper bounded by a composition of the leakage from all iterations.

For convex optimization with Lipschitz continuity, where the norms of gradients are uniformly bounded by some given constant, DP-SGD is known to produce an asymptotically tight privacy-utility tradeoff <cit.>.
However, it is, in general, impractical to assume Lipschitz continuity in tasks such as deep learning.
Either asymptotically or non-asymptotically, the study of practical implementations of DP-SGD with more realistic and specific assumptions remains very active <cit.> and demanding.
Much research effort has been dedicated to tackling the following two fundamental questions.
First, provided that we do not need to publish the intermediate computation results, how conservative is the privacy claim offered by DP-SGD?
Second, during practical implementation, how to properly select the training model and hyper-parameters?

Regarding the first question, many prior works <cit.> tried to empirically simulate what the adversary can infer from models trained by DP-SGD.
In particular, <cit.> examined the respective power of “black-box” and “white-box” adversaries, and suggested that a substantial gap between the DP-SGD privacy bound and the actual privacy guarantee  may exist.
Unfortunately, beyond DP-SGD, there are few known ways to produce, let along improve, rigorous DP analysis for a training process.
Most existing analyses need to assume either access to additional public data <cit.>, or strongly-convex loss functions to enable objective perturbation <cit.>.
Thus, for general applications, we still have to adopt the conservative DP-SGD analysis for the worst-case DP guarantees.

The second question is of particular interest to practitioners.
The implementation of DP-SGD is tricky as the performance of DP-SGD is highly sensitive to the selection of the training model and hyper-parameters.
The lack of theoretical analysis on gradient clipping makes it hard to find good parameters and optimize model architectures instructively, though many heuristic observations and optimizations on these choices are reported.
<cit.> showed how to find proper model architectures to balance learning capability and utility loss when the dataset is given.
[Most prior works report the best model and parameter selection by grid searching, where the private data is reused multiple times and the selection of parameters itself is actually sensitive.
This additional privacy leakage, partially determined by the prior knowledge on the training data is in general very hard to quantify, though in practice it might be small <cit.>.]
Recent work <cit.> also demonstrated empirical improvements through selecting the clipping threshold adaptively.
However, even with these efforts, there is still a long way to go if we want to practically train large neural networks with rigorous and usable privacy guarantees.
The biggest bottlenecks include

* the huge model dimension, which may be even larger than the size of the training dataset, while the magnitude of noise required for gradient perturbation is proportional to the square root of model dimension, and

* the long convergence time, which implies a massive composition of privacy loss and also forces DP-SGD to add formidable noise resulting in intolerable accuracy loss.

To this end, Tramer and Boneh in <cit.> argued that within the current framework, for most medium datasets ($<$ 500K datapoints) such as CIFAR10/100 and MNIST, the utility loss caused by DP-SGD will offset the powerful learning capacity offered by deep models.
Therefore, simple linear models usually outperform the modern Convolutional Neural Network (CNN) for these datasets.
How to privately train a model while still being able to enjoy the state-of-the-art success in modern machine learning is one of the key problems in the application of DP.

§.§ Our Strategy and Results

In this paper, we set out to provide a systematic study of DP-SGD from both theoretical and empirical perspectives to understand the two important but artificial operations: (1) iterative gradient perturbation and (2) gradient clipping.
We propose a generic technique, ModelMix, to significantly sharpen the utility-privacy tradeoff.
Our theoretical analysis also provides instruction on how to select the clipping parameter and quantify privacy amplification from other practical randomness.

We will stick to the worst-case DP guarantee without any relaxation, but view the private iterative optimization process from a different angle.
In most practical deep learning tasks, with proper use of randomness, we will have a good chance of finding some reasonable (local) minimum via SGD regardless of the initialization (starting point) and the subsampling <cit.>.
In particular, for convex optimization, we are guaranteed to approach the global optimum with a proper step size.
In other words, there are an infinite number of potential training trajectories [We will use training trajectory in the following to represent the sequences of intermediate updates produced by SGD.] pointing to some (local) minimum of good generalization, and we are free to use any one of them to find a good model.
Thus, even without DP perturbation, the training trajectory has potential entropy if we are allowed to do random selection.

From this standpoint, a slow convergence rate when training a large model is not always bad news for privacy.
This might seem counter-intuitive.
But, in general, slow convergence means that the intermediate updates wander around a relatively large domain for a longer time before entering a satisfactory neighborhood of (global/local) optimum. Training a larger model may produce a more fuzzy and complicated convergence process, which could compensate the larger privacy loss composition caused in DP-SGD.
We have to stress that our ultimate goal is to privately publish a good model, while DP-SGD with exposed updates is merely a tool to find a trajectory with analyzable privacy leakage.
The above observation inspires a way to find a better DP guarantee even under the conservative “whitebox” adversary model: can we utilize the potential entropy of the training trajectory while still bounding the sensitivity to produce rigorous DP guarantees?

To be specific, different from standard DP-SGD which randomizes a particular trajectory with noise, we aim to privately construct an envelope of training trajectories, spanned by the many trajectories converging to some (global/local) minimum, and randomly generate one trajectory to amplify privacy.
To achieve this, we must carefully consider the tradeoff between (1) controlling the worst-case sensitivity in the trajectory generalization and (2) the learning bias resultant from this approach.
We summarize our contributions as follows.

* We present a generic optimization framework, called ModelMix, which iteratively builds an envelope of training trajectories through post-processing historical updates, and randomly aggregates those model states before applying gradient descent.
We provide rigorous convergence and privacy analysis for ModelMix, which enables us to quantify  $(\epsilon, \delta)$-DP budget of our protocol.
The refined privacy analysis framework proposed can also be used to capture the privacy amplification of a large class of training-purpose-oriented operations commonly used in deep learning.
This class of operations include data augmentation <cit.> and stochastic gradient Langevin dynamics (SGLD) <cit.>, which cannot produce reasonable worst-case DP guarantees by themselves.

* We study the influence of gradient clipping in private optimization and present the first generic convergence rate analysis of clipped DP-SGD in deep learning.
To our best knowledge, this is the first analysis of non-convex optimization via clipped DP-SGD with only mild assumptions on the concentration of stochastic gradient.
We show that the key factor in clipped DP-SGD is the sampling noise [The noise corresponds to using a minibatch of samples to estimate the true full-batch gradient.] of the stochastic gradient.
We then demonstrate why implementation of DP-SGD by clipping individual sample gradients can be unstable and sensitive to the selection of hyper-parameters. Those analyses can be used to instruct how to select hyper-parameters and improve network architecture in deep learning with DP-SGD.

* ModelMix is a fundamental improvement to DP-SGD, which can be applied to almost all applications together with other advances in DP-SGD, such as low-rank or low-dimensional gradient embedding <cit.> and fine-tuning based transfer learning <cit.> (if additional public data is provided).
In our experiments, we focus on computer vision tasks, a canonical domain for private deep learning.
We evaluate our methods on CIFAR-10, FMNIST and SVHN datasets using various neural network models and compare with the state-of-the-art results.
Our approach improves the privacy/utility tradeoff significantly. For example, provided a privacy budget $(\epsilon=8,\delta=10^{-5})$, we are able to train Resnet20 on CIFAR10 with accuracy $70.4\%$ compared to $56.1\%$ when applying regular DP-SGD.
As for private transfer learning on CIFAR10, we can improve the $(\epsilon=2,\delta=10^{-5})$-DP guarantee in <cit.> to $(\epsilon=0.64,\delta=10^{-5})$ producing the same $92.7\%$ accuracy.

The remainder of this paper is organized as follows.
In Section <ref>, we introduce background on statistical learning, differential privacy and DP-SGD.
In Section <ref>, we formally present the ModelMix framework, whose utility in both convex and non-convex optimizations is studied in Theorem <ref> and Theorem <ref>,  respectively.
In Section <ref>, we show how to efficiently compute the amplified $(\epsilon, \delta)$ DP security parameters in Theorem <ref> and a non-asymptotic amplification analysis is given in Theorem <ref>.
Further experiments with detailed comparisons to the state-of-the-art works are included in Section <ref>.
Finally, we conclude and discuss future work in Section <ref>.

§.§ Related Works

Theoretical (Clipped) DP-SGD Analysis: When DP-SGD was first proposed <cit.>, and in most theoretical studies afterwards <cit.>, the objective loss function is assumed to be $L$-Lipschitz continuous, where the $L_2$ norm of the gradient is uniformly bounded by $L$.
This enables a straightforward privatization on SGD by simply perturbing the gradients.
In particular, for convex optimization on a dataset of $n$ samples, a training loss of $\Theta({\sqrt{d\log(1/\delta)}}/{(\epsilon n)})$ is known to be tight under an $(\epsilon, \delta)$ DP guarantee <cit.>.

However, it is hard to get a (tight) Lipschitz bound for general learning tasks.
A practical version of DP-SGD was then presented in <cit.>, where the Lipschitz assumption is replaced by gradient clipping to ensure bounded sensitivity.
This causes a disparity between the practice and the theory as classic results <cit.> assuming bounded gradients cannot be directly generalized to clipped DP-SGD.
Some existing works tried to narrow this gap by providing new analysis.
<cit.> presented a convergence analysis of smooth optimization with clipped SGD when the sampling noise in stochastic gradient is bounded.
But <cit.> requires the clipping threshold $c$ to be $\Omega(T)$ where $T$ is the total number of iterations.
This could be a strong requirement, as in practice, the iteration number $T$ can be much larger than the constant clipping threshold $c$ selected.
<cit.> relaxed the requirement with an assumption that the sampling noise is symmetric. <cit.> studied the special case where clipped DP-SGD is applied to generalized linear functions.
In this paper, we give the first generic analysis of clipped DP-SGD with only mild assumptions on the concentration property of stochastic gradients.

Assistance with Additional Public Data:
When additional unrestricted (unlabeled) public data is available, an alternative model-agnostic approach is Private Aggregation of Teacher Ensembles (PATE) <cit.>.
PATE builds a teacher-student framework, where private data is first split into multiple (usually hundreds) disjoint sets, and a teacher model is trained over each set separately.
Then, one can apply those teacher models to privately label public data via a private majority voting.
Those privately labeled samples are then used to train a student model, as a postprocessing of labeled samples.
Another line of works considers improving the noise added in DP-SGD with public data.
For example, in private transfer learning, we can first pretrain a large model with public data and then apply DP-SGD with private data to fine-tune a small fraction of the model parameters <cit.>.
However, both PATE and private transfer learning have to assume a large amount of public data. Another idea considers the projection of the private gradient into a low-rank/dimensional subspace, approximated by public samples, to reduce the magnitude of noise <cit.>.
When the public samples are limited, DP-SGD with low-rank gradient representation usually outperforms the former methods.

Except for PATE, our methods can be, in general, used to further enhance those state-of-the-art DP-SGD improvements with public data.
For example, using 2K ImageNet public samples, the low-rank embedding method in <cit.> can train Resnet20 on CIFAR10 with $79.1\%$ accuracy at a cost of $(\epsilon=111.2,\delta=10^{-5})$-DP; while ModelMix can improve the DP guarantee to $(\epsilon=6.1,\delta=10^{-5})$-DP with the same accuracy as shown in Section <ref>.

§ PRELIMINARIES

Empirical Risk Minimization: In statistical learning, the model to be trained is commonly represented by a parameterized function $f(w,x): (\mathcal{W},\mathcal{X}) \to \mathbb{R}$, mapping feature $x$ from input domain $\mathcal{X}$ into an output (prediction/ classification) domain. In the following, we will always use $d$ to represent the dimensionality of the parameter $w$, i.e., $w \in \mathbb{R}^d$.
For example, one may consider $f(w,x)$ as a neural network with a sequence of linear layers connected by non-linear activation layers, and $w$ represents the weights to be trained. Given a set $\mathcal{D}$ of $n$ samples $\{(x_i, y_i), i=1,2,...,n\}$, we define the problem of Empirical Risk Minimization (ERM) for some loss function $l(\cdot,\cdot)$ as follows,

\begin{equation}\small
\label{obg_loss_1}
\min_{w} F(w) = \min_{w} \frac{1}{n} \cdot \sum_{i=1}^n l(f(w,x_i),y_i).
\end{equation}

For convenience, we simply use $f(w, x_i, y_i)$ to denote the objective loss function $l(f(w,x_i),y_i)$ in the rest of the paper.
Below, we formally introduce the definitions of Lipschitz continuity, smoothness and convexity, which are commonly used in optimization research.

A function $g$ is $L$-Lipschitz if for all $w, w'\in\mathcal{W}$, $|g(w)-g(w')|\leq L\|w-w'\|_2$.

A function $g$ is $\beta$-smooth on $\mathcal{W}$ if for all $w, w'\in \mathcal{W}$, $g(w')\leq g(w)+\langle \nabla g(w), w'-w \rangle+\frac{\beta}{2}\|w'-w\|_2^2$.

A function $g$ is convex on $\mathcal{W}$ if for all $w, w'\in \mathcal{W}$ and $t \in (0,1)$, $f(t w + (1-t) w')\leq t g(w)+ (1-t) g(w')$.

In the following, we will simply use $\|\cdot\|$ to denote the $l_2$ norm unless specified otherwise.

Differential Privacy (DP): We first formally define $(\epsilon, \delta)$-DP and $(\alpha, \epsilon)$-Rényi DP as follows.

Given a data universe $\mathcal{X}^*$, we say that two datasets $\mathcal{D},\mathcal{D}'\subseteq \mathcal{X}^*$ are neighbors, denoted as $\mathcal{D} \sim \mathcal{D}'$, if $\mathcal{D} = \mathcal{D}' \cup s$ or $\mathcal{D}' = \mathcal{D} \cup s$ for some additional datapoint $s$. A randomized algorithm $\mathcal{A}$ is said to be $(\epsilon,\delta)$-differentially private (DP) if for any pair of neighboring datasets $\mathcal{D},\mathcal{D}'$ and any event $S$ in the output space of $\mathcal{A}$, it holds that

\begin{equation*}
\mathbb{P}(\mathcal{A}(\mathcal{D})\in S)\leq e^{\epsilon}\cdot \mathbb{P}(\mathcal{A}(\mathcal{D}')\in S)+\delta.
\end{equation*}

A randomized algorithm $\mathcal{A}$ satisfies $(\alpha, \epsilon)$-Rényi Differential Privacy (RDP), $\alpha >1$, if for any pair of neighboring datasets $\mathcal{D} \sim \mathcal{D}'$,

\[ \small
\epsilon \geq \mathsf{D}_{\alpha}(\mathcal{M}(\mathcal{D}) \| \mathcal{M}(\mathcal{D}') ).
\]

\begin{equation}\small
\label{Rényi-diver}
\mathsf{D}_{\alpha}(\mathsf{P} \| \mathsf{Q}) = \frac{1}{\alpha-1} \log \int \mathsf{q}(o) (\frac{\mathsf{p}(o)}{\mathsf{q}(o)})^{\alpha} ~d o,
\end{equation}

represents $\alpha$-Rényi Divergence between two distributions $\mathsf{P}$ and $\mathsf{Q}$ whose density functions are $\mathsf{p}$ and $\mathsf{q}$, respectively.

In Definition <ref> and <ref>, if two neighboring datasets $\mathcal{D}$ and $\mathcal{D}'$ are defined in a form that $\mathcal{D}$ can be obtained by arbitrarily replacing an datapoint in $\mathcal{D}'$, then
they become the definitions of bounded DP <cit.> and RDP, respectivaly. In this paper, we adopt the unbounded DP version to match existing DP deep learning works <cit.> with a fair comparison.

In practice, to achieve meaningful privacy guarantees, $\epsilon$ is usually selected as some small one-digit constant and $\delta$ is asymptotically $O(1/|\mathcal{D}|) = O(1/n)$.
To randomize an algorithm, the most common approaches in DP are Gaussian or Laplace Mechanisms <cit.>, where a Gaussian or Laplace noise proportional to the sensitivity is added to perturb the algorithm's output.
In many applications, including the DP-SGD analysis, we need to quantify the cumulative privacy loss across sequential queries of some differentially private mechanism on one dataset. The following theorem provides an upper bound on the overall privacy leakage.

For any $\epsilon>0$ and $\delta \in (0,1)$, the class of $(\epsilon, \delta)$-differentially private mechanisms satisfies $(\tilde{\epsilon}, T\delta + \tilde{\delta})$-differential privacy under $T$-fold adaptive composition for any $\tilde{\epsilon}$ and $\tilde{\delta}$ such that

\[
\tilde{\epsilon} = \sqrt{2T\log(1/\tilde{\delta})}\cdot\epsilon + T\epsilon(e^{\epsilon}-1).
\]

For any $\alpha>1$ and $\epsilon>0$, the class of $(\alpha, \epsilon)$-RDP mechanisms satisfies $(\tilde{\epsilon}, \tilde{\delta})$-differential privacy under $T$-fold adaptive composition for any $\tilde{\epsilon}$ and $\tilde{\delta}$ such that

\[
\tilde{\epsilon} = T\epsilon - \log(\tilde{\delta})/(\alpha-1).
\]

Theorem <ref> provides a good characterization on how the privacy loss increases with composition.
For small $(\epsilon,\delta)$, we still have an $\tilde{O}(\sqrt{T}\epsilon, T\delta)$ DP guarantee after a $T$ composition.
In practice using RDP, Theorem <ref> usually produces tighter constants in the privacy bound.

DP-SGD: (Stochastic) Gradient Descent ((S)GD) is a very popular approach to optimize a function.
Suppose we try to solve the ERM problem and minimize some function $F(w) = \frac{1}{n} \sum_{i=1}^n f(w,x_i, y_i)$.
SGD can be described as the following iterative protocol.
In the $(k+1)$-th iteration, we apply Poisson sampling, i.e., each datapoint is i.i.d. sampled by a constant rate $q$, and a minibatch of $B_{k}$ samples is produced from the dataset $\mathcal{D}$, denoted as $S_k$.
We calculate the stochastic gradient as

\begin{equation}
\label{stochastic gradient}
G_{k} \gets \sum_{(x_i, y_i) \in S_k} \nabla f(w_{k},x_i,y_i).
\end{equation}

Then, a gradient descent update is applied using

\begin{equation}
\label{SGD}
w_{k+1} = w_{k} - \eta \cdot G_k,
\end{equation}

for some stepsize $\eta$.
In particular, if the minibatch is selected to be the full batch, i.e., $S_k = \mathcal{D}$, then Equation (<ref>) becomes the standard gradient descent procedure.

We make the following assumption regarding the sampling noise $\| \nabla f(w,x,y) - \nabla F(w)\|$ when the minibatch size equals $1$.
This assumption, which will be shown to be necessary in Example <ref>,  will be used in Theorem <ref> when we derive the concrete convergence rate for clipped SGD.
[Stochastic Gradient of Sub-exponential Tail]

There exists some constant $\kappa>0$ such that for any $w$, if we randomly select a datapoint $(x,y)$ from $\mathcal{D}$, then

\[\Pr(\| \nabla f(w,x,y) - \nabla F(w)\| \geq t) \leq e^{-t/\kappa}.\]

In Assumption <ref>, a larger $\kappa$ implies stronger concentration, i.e., a faster decaying tail of the stochastic gradient.
The modification from GD/SGD to its corresponding DP version is straightforward.
When the loss function $f$ is assumed to be $L$-Lipschitz <cit.>, i.e., $\|\nabla f(w, x_i,y_i)\| \leq L$ for any $w$, the worst-case sensitivity in Equation (<ref>) is bounded by $\eta L$ in each iteration.
One can derive a tighter bound <cit.> using existing results on the privacy amplification from sampling <cit.>.
Thus, SGD can be made private via iterative perturbation by replacing Equation (<ref>) with the following:

\begin{equation}
\label{dpsgd_without_clipping}
w_{k+1} = w_{k} - \eta\cdot (G_k + \Delta_{k+1}),
\end{equation}

where $\Delta_{k}$ is the noise for the $k$-th iteration.
For example, if we want to use the Gaussian Mechanism to ensure $(\epsilon, \delta)$-DP when running $T$ iterations, then $\Delta_{k+1}$ can be selected to be i.i.d. generated from

\[
\Delta_{k+1} \leftarrow \mathcal{N}(\bm{0}, O({L^2T\log(1/\delta)\over \epsilon^2}) \cdot \bm{I}_d).
\]

Here, $\bm{I}_d$ represents the $d \times d$ identity matrix.

However, when we do not have the Lipschitz assumption, an alternative is to force a limited sensitivity through gradient clipping.
Following the same notations as before, we describe DP-SGD with per-sample gradient clipping <cit.> as follows,

\begin{equation}
\begin{aligned}
\label{dpsgd_with_cliping}
& G_k \gets \sum_{(x_i, y_i) \in S_k} \mathsf{CP}\big(\nabla f(w_k,x_i,y_i), c\big);\\
& w_{k+1} = w_{k} - \eta\cdot (G_k + \Delta_{k+1}).
\end{aligned}
\end{equation}

Here, $\mathsf{CP}(\cdot, c)$ represents a clipping function of threshold $c$,

$$\mathsf{CP}(\nabla f(w,x,y), c) = \nabla f(w,x,y) \cdot \min\{1, {c \over \|\nabla f(w,x,y) \|} \}.$$

With clipping, the $l_2$ norm of each per-sample gradient is bounded by $c$.
Thus, the clipping threshold $c$ virtually plays the role of the Lipschitz constant $L$ in clipped SGD for privacy analysis.

§ MODELMIX

In this section, we formally introduce ModelMix, and explain how it sharpens the utility-privacy trade-off in DP-SGD.

§.§ Intuition

We begin with the following observation.
Suppose we run SGD twice on a least square regression $F(w) = 1/n \cdot \sum_{i=1}^n \|\langle w, x_i\rangle -y_i\|^2$ for $T$ iterations and obtain two training trajectories $\bm{w}=(w_1,w_2,...,w_T)$ and $\bm{w}'=(w'_1,w'_2,...,w'_T)$.
Suppose both $w_T$ and $w'_T$ are $\sigma$-close to the optimum $w^* = \arg \min_{w} F(w)$, i.e.,

\[
\|w_T-w^*\|\leq \sigma \text{~and~} \|w'_T-w^*\|\leq \sigma.
\]

Due to the linearity of gradients in least square regression, if we mix $\bm{w}$ and $\bm{w}'$ to get $w''_k = \alpha w_{k} + (1-\alpha) w'_{k}$ ($k=1,2,...,T$) for some weight $\alpha \in  (0,1)$, then we produce a new SGD trajectory $\bm{w}''$ where $w''_T$ is also $\sigma$-close to $w^*$.

This simple example gives us two inspirations.
First, as mentioned earlier, with different randomness in initialization and subsampling, the training trajectory to find an optimum is not unique.
Even in DP-SGD where we need to virtually publish the trajectory for analytical purpose, we are not restricted to expose a particular one.
Second, and more important, this means that we have more freedom to randomize the SGD process.
In the $k$-th iteration, instead of following the regular SGD rule in Equation (<ref>) where we simply start from the previously updated state $w_{k-1}$, we can randomly mix $w_{k-1}$ with some other $w'_{k-1}$ from another reasonable training trajectory, and then move to a new trajectory to proceed.

However, the above idea to randomly mix the trajectories cannot be directly implemented as we do not have any prior knowledge on what good trajectories look like.
Any training trajectory generated by the private dataset is sensitive and potentially creates privacy leakage.
Thus, we need to be careful about how we generate the needed envelope.
Recall that we use envelope to describe the space spanned by the mixtures of training trajectories. Provided the property that DP is immune to post-processing, we consider approximating the trajectory envelope using intermediate states already published. This allows us to privately construct the envelope and advance the optimization, simultaneously.

§.§ Algorithm and Observations

We start with a straw-man solution where we virtually run DP-SGD to alternately train two models in turns.
We initialize two states $\tilde{w}^1_0$ and $\tilde{w}^2_0$ with respect to (w.r.t) the parameterized function we aim to optimize.

At any odd iteration, i.e., iteration $2k + 1$ ($k \ge 0$), we randomly generate $\bm{\alpha}_{2k+1} \in (0,1)^d$ whose coordinate is i.i.d. uniform in $(0,1)$.
We then update the state of the first model $\tilde{w}^1$ as

\begin{equation}
\label{distributed-manner-1}
\bar{w}^1_{k} \gets
\bm{\alpha}_{2k+1} \circ \bar{w}^1_{k-1}  + (\bm{1}-\bm{\alpha}_{2k + 1}) \circ \bar{w}^2_{k-1} - \eta \big( \nabla F(\tilde{w}^1_{k-1})+\Delta_{2k+1} \big),
\end{equation}

where $\circ$ represents the Hadamard product and $\Delta$ represents the noise added.
Essentially, we mix the two states $\bar{w}^1_{k-1}$ and $\bar{w}^2_{k-1}$ coordinate-wise.
Similarly, at any even iteration, i.e., iteration $2(k+1)$, we randomly generate some $\bm{\alpha}_{2k+2}$ and update the state of the second model $\tilde{w}^2$ as

\begin{equation}
\label{distributed-manner-2}
\bar{w}^2_{k} \gets \bm{\alpha}_{2k+2} \circ \bar{w}^1_{k}  + (\bm{1}-\bm{\alpha}_{2k + 2}) \circ \bar{w}^2_{k-1} - \eta \big( \nabla F(\tilde{w}^2_{k-1})+\Delta_{2k+2} \big).
\end{equation}

We take turns mixing and updating two training trajectories privately using (<ref>) and (<ref>).
At the high level, this is similar to a distributed GD, where two agents collaboratively train the model.

However, in a centralized scenario, it is unnecessary to artificially create two models and  train them in a distributed manner, meaning that the efficiency of the above mixture method is not optimal.
In the following, we propose an improved method where in the $k$-th iteration, we simply post-process on $w_{k-1}$ and $w_{k-2}$, the updates already privately generated from the last two iterations, instead of two virtual models to approximate an envelope of training trajectory.
The formal description of ModelMix is shown in Algorithm <ref>. We provide an illustration of how Algorithm <ref> works in Fig. <ref>.

Differentially Private Stochastic Gradient Descent with ModelMix
[1]
Input: Objective function $F(w)= {1\over n} \cdot \sum_{i=1}^n f(w, x_i, y_i)$, dataset $\mathcal{D}=\{(x_i,y_i), i=1,2,...,n\}$, sampling rate $q$, step size $\eta$, clipping threshold $c$, total number of iterations $T$, mixing thresholds $\tau_{[1:T]}$, initialized model states $w_{-1}, w_{0} \in \mathbb{R}^d$ and noise sequence $\Delta_{[1:T]}$.
Through i.i.d. sampling with a constant rate $q$, produce a minibatch $S_{k}$ of $B_k$ samples from $\mathcal{D}$ and calculate the stochastic gradient
\[G_{k-1} =  \sum_{(x_i,y_i) \in S_k} \mathsf{CP}\big(\nabla f(w_{k-1},x_i,y_i), c \big).\]

$\bm{\alpha}_k(j)\leftarrow \mathcal{U}[0,1]$, a uniform distribution in [0,1].
$|{w}_{k-1}(j) - {w}_{k-2}(j)| < \tau_k $
${w}_{k-1}(j) \gets {w}_{k-1}(j)+ \text{sign}({w}_{k-1}(j)-{w}_{k-2}(j)) \cdot \tau_k/2$;
${w}_{k-2}(j) \gets {w}_{k-2}(j)- \text{sign}({w}_{k-1}(j)-{w}_{k-2}(j)) \cdot \tau_k/2$.
Update the weight as follows:

\begin{equation}
\label{main_update}
\begin{aligned}
\bm{\alpha}_k(j)\cdot{w}_{k-1}(j)  + (1-\bm{\alpha}_k(j))\cdot {w}_{k-2}(j) \\
& -\eta\cdot (G_{k-1}(j)+\Delta_k(j)).
\end{aligned}
\end{equation}

Output: $w_{T}$.

An illustration of how ModelMix works.
Each node represents an intermediate state.
The dashed lines represent the mixing and the arrows represent gradient descent. The green area represents the approximated trajectory envelope.

From a more concise optimization standpoint, the mixed state in expectation is the average of the last two iterations' states $(w_{k-1}+w_{k-2})/2$.
Therefore, in the expectation of Equation (<ref>), we are essentially optimizing the original objective function $F(w)$ plus a proximal term in a form

\begin{equation}
\label{proximal}
F(w) + \|w - w_{k-2} \|^2/{4},
\end{equation}

whose gradient at $w_{k-1}$ is $\nabla F(w_{k-1}) + (w_{k-1}-w_{k-2})/2$.
Thus, ModelMix can also be viewed as introducing a randomized proximal term into the original objective function $F(w)$. ModelMix operation $\bm{\alpha}_{k}w_{k-1}+ (\bm{1}-\bm{\alpha}_{k-1})w_{k-2}$ averages out the noise added in the previous iterations, which makes the convergence more stable. This enables us to apply a larger step size during training, as shown later in Fig.  <ref>.

As will be shown later in Theorem <ref>, the privacy amplification is determined by the distance between $w_{k-1}$ and $w_{k-2}$.
Therefore, for a privacy analysis purpose, in steps 5-8 of Algorithm <ref>, we artificially ensure that the coordinate-wise distance between the states $w_{k-1}$ and $w_{k-2}$ is at least some parameter $\tau_k$.
In practice, the gap between the two states $w_{k-1}$ and $w_{k-2}$ already exists even without these operations, meaning that ModelMix naturally enjoys a privacy amplification.
We only add those steps to enforce a worst-case lower bound on the coordinate-wise distance so that we can quantify the privacy amplification in a clean way.

A natural question that follows is how do the parameters $\tau_{[{1:T]}}$ across $T$ iterations affect privacy and utility?
$\tau_k$ captures the distance between the two training trajectories, or in other words, the volume of the approximated envelope.
Ideally, we would like $\tau_{[{1:T]}}$ to match the size of the true envelope so that the impact of those artificial steps is minimized.
Empirically, a harder learning task often has a slower convergence rate, and thus also has an envelope of larger volume.
This allows us to select a larger $\tau_k$ and provide stronger privacy amplification.
On the other hand, if the learning task is simple and we overestimate the size of the envelope with a large $\tau_k$, then the learning rate might be compromised.
This analysis is supported by Fig. <ref>, where we test the efficiency of ModelMix in two different tasks with various setups.

In Fig. <ref> (a) and (b), we implement non-private SGD with or without ModelMix to train Resnet20 <cit.> on CIFAR10 [https://www.cs.toronto.edu/ kriz/cifar.html] and SVHN [http://ufldl.stanford.edu/housenumbers/], the two benchmark datasets.
We set the sampling rate $q$ to be $0.05$ and run for 2,000 iterations.
Compared to CIFAR10, classification on SVHN is an easier task.
We fix the mixing threshold $\tau_k = \tau$ to be proportional to the step size $\eta$ in gradient descent (Equation (<ref>)), specifically, $0.1\eta$ and $0.05\eta$.

In both Fig. <ref> (a) and (b), different setups perform similarly in the initial stage (the first 400 iterations).
However, in later epochs, as the training trajectories approach the (local) optimum, setups with larger $\tau$ suffer heavier losses.
This phenomenon is clearer in Fig. <ref> (b) which trains on the easier SVHN task.
This matches our previous analysis: the size of the trajectory envelope gets smaller when (1) the learning task is easier or (2) we approach the optimum.
On the other hand, with proper selection of $\tau=0.05\eta$, in the non-private case ModelMix comes with almost no additional utility loss.

In Fig. <ref> (c) and (d), we implement the private case where we will see how ModelMix strengthens the robustness of DP-SGD.
We clip the per-sample gradient norm down to $8$ and add the same amount of Gaussian noise to both DP-SGD cases, with or without ModelMix.
This ensures an $(\epsilon=8, \delta=10^{-5})$ guarantee for regular DP-SGD.
As shown later in Theorem <ref>, ModelMix can achieve a much better DP guarantee under the same setup.
But here, we focus on the performance of ModelMix given the same noise as regular DP-SGD to provide a clear comparison.
With the virtual proximal term in Equation (<ref>), ModelMix allows us to use a larger stepsize resulting in a faster learning rate.
In  Fig. <ref> (c) and (d), we set the stepsize $\eta=0.2$ for ModelMix, and compare to the regular DP-SGD with $\eta=0.1$ (black line) and $\eta=0.2$ (red line). The larger stepsize $\eta=0.2$ worsens the performance of regular DP-SGD, whereas with proper selection of $\tau=0.05\eta$ (green line), ModelMix allows the application of larger stepsize and slightly out-performs regular DP-SGD.

Performance comparison between DP-SGD with/without ModelMix (MM) provided with the same noise.

§.§ Utility Guarantee in (non)-convex Optimization and Effect of Gradient Clipping

In this subsection, we will provide formal utility analysis of Algorithm <ref>.
As a warm-up, in Theorem <ref>, we start with convex optimization with a Lipschitz and smooth assumption, and show Algorithm <ref> enjoys an $1/\sqrt{T}$ convergence rate.
Specifically, after $T$ iterations, if we set the average

\[
\bar{w} = \frac{\sum_{k=1}^T (w_{k-1}+w_{k-2})}{2T},
\]

where $w_{k}$ represents the weight in the $k$-th iteration, then the utility loss $F(\bar{w}) - F(w^*)$ is bounded by $O(1 / \sqrt{T})$. Based on Theorem <ref>, we then move to the non-convex case without a Lipschitz assumption in Theorem <ref>, and show a generic convergence rate.

Suppose $f(w,x,y)$ is $L$-Lipschitz and $\beta$-smooth convex.
We set the clipping threshold $c = L$.
For any $\gamma>0$, if we set $\eta = \gamma /(nq\sqrt{T})$, then Algorithm <ref> satisfies

\begin{equation}
\begin{aligned}
\mathbb{E} \big[ F(\bar{w}) - F(w^*) \big] \\
& \leq \frac{3\mathcal{W}^2_0+ \sum_{k=1}^Td \tau^2_{i}/12}{2\gamma \sqrt{T}} + \frac{\gamma(L^2/q^2+\mathbb{E}[\|\Delta\|^2]/(n^2q^2))}{\sqrt{T}} + \beta \cdot \\
& \big( \frac{12\mathcal{W}^2_0}{8T} + \frac{2\gamma \mathcal{W}_0(L+\mathbb{E}[\|\Delta\|]/n)}{q T^{3/2}} + \frac{11\gamma^2(L^2+\mathbb{E}[\|\Delta\|^2]/n^2)}{8Tq^2} \big) \\
& = O(\frac{\mathcal{W}^2_0+ \sum_{k=1}^Td\tau^2_{k}}{\gamma \sqrt{T}} + \frac{\gamma(L^2+\mathbb{E}[\|\Delta\|^2]/n^2)}{q^2\sqrt{T}})
\end{aligned}
\label{convex_convergence}
\end{equation}

where  $\mathcal{W}_0 = \sup_{w}\|w-w^*\|$ denotes the initial divergence and $\mathbb{E}[\| \Delta\|^2] = \mathbb{E}[\| \Delta_{k}\|^2]$ denotes the variance of noise added in each iteration.

See Appendix <ref>.

From Equation (<ref>), one can see that with ModelMix, Algorithm <ref> still enjoys an $O({1}/{\sqrt{T}})$ convergence rate when $\tau_k = O(\eta) = O(1/\sqrt{T})$.
Using Theorem <ref>, we can upper bound the utility loss of Algorithm <ref>. Since ModelMix can be seen as post-processing operations on DP-SGD, its privacy guarantee is at least as good as the privacy guarantee of DP-SGD.
If we set $\Delta_{k}$ to be Gaussian noise generated from $\mathcal{N}(0,O(q^2L^2T\log(1/\delta)/\epsilon^2) \cdot \bm{I}_{d} )$ <cit.>, then by Theorem <ref>, Algorithm <ref> satisfies an $(\epsilon, \delta)$-DP guarantee and its utility loss is upper bounded by $O(\sqrt{d}\log(1/\delta)L/(n\epsilon))$.
This asymptotically matches the classic results in <cit.>. However, we will present a more fine-tuned analysis in Section <ref> to show the randomness in ModelMix can significantly sharpen the composite privacy loss.

In the following, we no longer assume the objective loss function to be either convex or Lipschitz continuous, and we try to understand the effect of gradient clipping.
The following theorem gives a generic convergence analysis of ModelMix in non-convex optimization with clipped gradient.

Suppose the objective loss function $F(w)$ is $\beta$-smooth and satisfies Assumption <ref>, then there exists some constant $\psi > 0$ such that when the clipping threshold $c$ satisfies

\begin{equation}
\label{selection-of-c}
c \geq \max \{4\kappa\log(10), -\psi\kappa\log(\kappa)\log(\frac{\sqrt{d \log(1/\delta)}}{n\epsilon})\},
\end{equation}

then the convergence rate of Algorithm <ref>, which applies per-sample gradient clipping up to $c$ and enjoys an $(\epsilon, \delta)$-DP guarantee, satisfies

\begin{equation}
\begin{aligned}
\label{batch_main_1}
&  ~~~~~\mathbb{E} \Big[ \frac{\sum_{k=0}^{T-1}  \min\big\{ 9/20 \cdot  \|\nabla F(w_k)\|^2, c/20 \cdot \| \nabla F(w_k)\| \big \}}{T}\Big]  \\
& \leq (\frac{v}{2} +\frac{5}{2})\frac{c\sqrt{\mathcal{R}_{F}\frac{101}{12}\beta d \log(1/\delta)}}{n\epsilon} + \frac{28c\beta{d\log(1/\delta)}\tilde{\mathcal{W}}_0 }{12q (n\epsilon)^2} \\
&  ~~~~+  \frac{c d\log(1/\delta)\sqrt{\frac{101}{12}}\beta^{3/2}}{qn\epsilon\sqrt{\mathcal{R}_{F}}} \big(\frac{\sum_{k=1}^T d\tau^2_{k}}{12} + \frac{21\tilde{\mathcal{W}}_0^2}{24}\big),
\end{aligned}
\end{equation}

where $\mathcal{R}_{F} = \sup_{w}F(w)- \inf_{w} F(w)$, $v$ is some constant determined by the noise mechanism and $\tilde{\mathcal{W}}_0 = \|w_0-w_{-1}\|$ is the initial divergence. When $\tau_{k} = O(\eta)$, the right hand of  (<ref>) is $ O(\frac{c\sqrt{d \log(1/\delta)}}{n\epsilon})$.

See Appendix <ref>.

In Theorem <ref>, we do not assume the objective function $F(w)$ to be convex or Lipschitz, but only smooth with concentrated stochastic gradients  (Assumption <ref>).[In the proof of Theorem <ref>, essentially we only need a high-probability bound of the sampling noise in the stochastic gradient, which is also necessary (see Example <ref>). Therefore, one can relax Assumption <ref> of a subexponential tail to any proper concentration assumption that allows the derivation of a high probability bound of the sampling noise.
The utility loss is measured by the norm of the gradient $\|\nabla F(w_k)\|$ in Equation (<ref>), commonly considered in non-convex optimization <cit.>.
There are several interesting observations from Theorem <ref>.
First, in clipped DP-SGD, we should not simply consider the clipping threshold $c$ as a virtual Lipschitz constant.
When the gradient norm $\|\nabla F(w_k)\|$ is large, Equation (<ref>) suggests that
$${\sum_{k=0}^{T-1}\|\nabla F(w_k)\|} = O(T \cdot \frac{\sqrt{\log(1/\delta)d}}{n\epsilon}),$$

which is independent of the clipping threshold $c$ as long as $c$ satisfies Equation (<ref>).
Therefore, when the objective function is hard to optimize, or when the gradient is large during the initial epochs, we should select a large clipping threshold to minimize the effect of clipping and let the state $w_k$ approach some neighborhood domain of a (local) minimum fast.
When $\|\nabla F(w_k)\|$ becomes small, $c$ will reappear in the utility bound.
In this case, the equation becomes similar to the classic DP-SGD privacy-utility tradeoff <cit.> but replacing the Lipschitz constant $L$ by $c$, where

$$\frac{\sum_{k=0}^{T-1}\|\nabla F(w_k)\|}{T} = O(\frac{c\sqrt{\log(1/\delta)d}}{n\epsilon}).$$

Therefore, we want to select some properly small $c$.

The second and also more important issue worth mentioning is the restriction on $c$ in Equation (<ref>).
We require the norm of the sampling noise in the stochastic gradient to be smaller than the clipping threshold $c$ with sufficient probability.
Indeed, this is actually a necessary condition for clipped DP-SGD to work.
In the following example, we show that even without any noise perturbation, clipped SGD could fail when $c$ is much smaller than the magnitude of sampling noise.

Suppose the loss function is $f(w, x) = (w - x)^2 / 2$ and we have three samples $x_1=-20,x_2=-10$ and $x_3=90$.
The overall loss function is thus $F(w)= 1/3 \cdot \sum_{i=1}^3 (w-x_i)^2/2$.
Thus, $\nabla F(w)=w-20$ and the minimum is achieved when $w=20$.
Moreover, the standard deviation  of the per-sample stochastic gradient, where we randomly sample $x'\leftarrow \{x_1, x_2, x_3\}$ and compute the gradient of $f(w, x')$, is $49.7$ for any $w$.

We will show that if the clipping threshold is relatively small, for example, $c=1$, then the clipped gradient is dramatically different from the true gradient $\nabla F(w)$.
We first consider the case when the magnitude of true gradient is small, say at the optimum $w=20$ where $\nabla F(20)=0$.
At $w=20$, the expectation of the clipped stochastic gradient is $1/3$ rather than $0$.
On the other hand, when the magnitude of the true gradient is large, say $w=0$ and $\nabla F(0)=-20$, the expectation of the clipped stochastic gradient is still $1/3$.
But now it is in the opposite direction of the true gradient.

From Example <ref>, we know that if $c$ is not properly selected, clipped SGD could fail to converge regardless of the magnitude of the full-batch gradient. As a summary, Theorem <ref> and Example <ref> suggest the following two key observations on clipped DP-SGD:

* On hyper-parameter selection, to ensure stable convergence in general, we need the clipping threshold $c$ to be sufficiently larger than the sampling noise of stochastic gradient.

* On the other hand, to improve the performance of clipped DP-SGD, one promising direction is to reduce the sampling noise via network architecture optimization and data normalization. Indeed, many approaches supporting this goal are proposed and studied in deep learning research, for example, batchnorm layer <cit.>, where the gradient is evaluated and normalized by a group of samples. We defer a systematical study on improvement via gradient variance reduction to our future work.

§ PRIVACY GUARANTEES AND AMPLIFICATIONS

§.§ Numerical Calculation and Asymptotic Amplification

In this section, we first show how to numerically compute the $(\epsilon, \delta)$ bound via RDP and asymptotically analyze the amplification of ModelMix.
We use Gaussian Mechanism to produce the noise $\Delta_k \sim \mathcal{N}(0 , \sigma^2 \cdot \bm{I}_{d})$ across iterations.
The following Theorem <ref>  shows an efficiently-calculable upper bound of $(\epsilon, \delta)$ by measuring an $\alpha$-Rényi divergence between two one-dimensional distributions.
For simplicity, we fix $\tau_k = \tau$ in the following.

Suppose $\Delta_k$ is Gaussian noise generated from $\mathcal{N}(0 , \sigma^2 \cdot \bm{I}_{d})$.
We define two distributions $\mathsf{P}_0$ and $\mathsf{P}_1$,

\[
\mathsf{P}_0 = \mathcal{N}(0,\sigma)*\mathcal{U}[-\tau/(2\eta),\tau/(2\eta)],
\]

\[
\mathsf{P}_1 = \mathcal{N}(c,\sigma)*\mathcal{U}[-\tau/(2\eta),\tau/(2\eta)],
\]

where $*$ represents the convolution of two probability distributions.
In other words, sampling from $\mathsf{P}_0$ is equivalent to sampling from Gaussian distribution $\mathcal{N}(0,\sigma)$ and uniform distribution $\mathcal{U}[-\tau/(2\eta),\tau/(2\eta)]$ independently, then summing the results together.
Similarly, $\mathsf{P}_1$ is the convolution of $\mathcal{N}(c,\sigma)$ and $\mathcal{U}[-\tau/(2\eta),\tau/(2\eta)]$.

Algorithm <ref> for $T$ iterations with sampling rate $q$ satisfies $(\epsilon, \delta)$ DP for any $\epsilon$ and $\delta$ such that

$$\epsilon \leq \min_{\alpha \in \mathbb{Z}, \alpha>1} T\mathsf{D}_{\alpha}\big( (1-q)\mathsf{P}_0 + q \mathsf{P}_1 \| \mathsf{P}_0 \big) + {\log(1/\delta)\over (\alpha-1)}.$$

See Appendix <ref>.

Let $\mathcal{D}$ and $\mathcal{D}'$ be two neighboring datasets, where $(\tilde{x}, \tilde{y}) = \mathcal{D}-\mathcal{D}'$ is the differing element.
Theorem <ref> shows that the worst output divergence between $\mathcal{D}$ and $\mathcal{D}'$ occurs when the gradient $\nabla f(w_k,\tilde{x}, \tilde{y})$ is a vector whose Hamming weight is 1, such as $(c,0,...,0)$.
In other words, given an $l_2$-norm budget of $c$, the worst-case $\mathsf{D}_{\alpha}$ divergence happens when the gradient of $(\tilde{x}, \tilde{y})$ is concentrated on a single dimension. This is different from standard (subsampled) Gaussian Mechanism, where the randomization is only due to the Gaussian noise and subsampling <cit.>.
We will exploit this property to design more fine-tuned gradient clipping for ModelMix in Section <ref>. In the following, we give an asymptotic analysis on the privacy amplification of ModelMix.

If the sampling rate $q$ is some constant, the noise $\Delta_{k}$ is selected to be Gaussian noise $\mathcal{N}(\bm{0},\sigma^2 \bm{I}_d)$ and $\tau$ is sufficiently large, then Algorithm
<ref> with $T$ iterations satisfies $(\epsilon, \delta)$-DP guarantee for any $\epsilon$ and $\delta$ such that

\begin{equation}
\begin{aligned}
\epsilon =  \tilde{O}\big(\frac{\eta T(c+\sigma)}{{\tau}n} & + \sqrt{
\frac{\eta T\log(1/\tilde{\delta})}{\tau}\cdot (\frac{c^4}{n^4\sigma^3}+ \frac{c^2}{n^2\sigma} + \sigma)  }\big),\\
\delta & = \tilde{O}(T\cdot \frac{(L/n)+\sigma}{\tau} + \tilde{\delta}). \footnote{\hl{In the case of clipped gradient, one can simply replace $L$ with $c$},}
\end{aligned}
\end{equation}

where $\tilde{\delta}$ can be any value within $(0,1)$.

In addition, if in Algorithm <ref>, the per-sample gradient is clipped up to $c$ in $l_1$ norm and the noise $\Delta_{k}$ is selected to be Laplace noise, then Algorithm <ref> satisfies $(\epsilon, \delta)$-DP, where

\begin{equation}
\epsilon =   \tilde{O}\big(\frac{\eta T\epsilon_0(e^{\epsilon_0}-1)}{\tau} + \epsilon_0 \sqrt{\frac{\eta T \log(1/\delta)}{\tau}}\big),
\end{equation}

where $\epsilon_0$ represents the $\epsilon$ privacy loss of a single iteration of (<ref>) without ModelMix (set $\bm{\alpha}_k$ to be constant).

See Appendix <ref>.

We can compare between Theorem <ref> and Theorem <ref> (the case of regular DP-SGD without ModelMix).  If a single iteration of a DP-SGD is $\epsilon_0$-differentially-private, then
classic advanced composition (Theorem <ref>) suggests that the composition of $T$ iterations produces $O(\epsilon_0\sqrt{T\log(1/\delta)},\delta)$-DP when $\epsilon_0$ is small.
Theorem <ref> states that under the same setup, if the DP-SGD is further incorporated with ModelMix, then this composition becomes $\tilde{O}(\epsilon_0\sqrt{T\log(1/\delta)/\tau}, \delta)$, where the $\epsilon$ term decreases by a factor of $\tilde{O}(1/\sqrt{\tau})$.

§.§ Clipping with $l_2$ and $l_{\infty}$ Sensitivity Guarantee

$l_2$ and $l_{\infty}$ Norm Gradient Clipping
[1]
Input: Individual gradient $\nabla f(w, x, y)$, $l_2$ norm clipping threshold $c$, $l_{\infty}$ norm truncation parameter $p \in \mathbb{Z}^{+}$.
Clip $\nabla f(w,x,y) \gets \mathsf{CP}(\nabla f(w,x,y), c)$
$j=1,2, ... , d$
$|\nabla f(w, x, y) (j)| > \frac{c}{\sqrt{p}}$
$\nabla f(w, x, y) (j) \gets
\text{sign} (\nabla f(w, x, y) (j)) \frac{c}{\sqrt{p}}$
Output: $\nabla f(w, x, y).$

In this subsection, we show that in the framework of ModelMix with Gaussian Mechanism, one can further strengthen the privacy amplification if both the $l_2$ and $l_{\infty}$ norm sensitivity can be guaranteed.
We first present a gradient clipping algorithm, described in Algorithm <ref>.
Algorithm <ref> is a simple generalization of standard $l_2$-norm clipping operator $\mathsf{CP}(\cdot, c)$. For an individual gradient $\nabla f(w,x,y)$,  we first clip the gradient up to $c$ in $l_2$ norm. After that, we truncate each coordinate such that its absolute value does not exceed $c/\sqrt{p}$ for some constant integer $p \geq 1$. As a consequence, Algorithm <ref> ensures that the $l_2$ and $l_{\infty}$ norm of clipped gradient is upper bounded by $c$ and $c/\sqrt{p}$, respectively.
From our empirical observations, compared to $l_2$-norm clipping on $c$, DP-SGD is much less sensitive to $l_\infty$-norm truncation on $p$. The reason behind this is that in practice the gradients obtained are rarely concentrated on few coordinates and thus $l_\infty$ truncation with large $p$, even in the hundreds, will hardly change the geometry of gradients. Roughly speaking, $p$ captures the number of significant coordinates in the gradient, where the  gradient dimension could be millions in deep learning. However, such truncation enables us to derive a stronger amplification bound with ModelMix, as shown in  Corollary <ref>.

Under the same setup of Theorem <ref>, if we further ensure that the $l_{\infty}$ norm of each individual gradient in Algorithm <ref> is bounded by $c/\sqrt{p}$ for some $p \in \mathbb{Z}^+$, then it satisfies $(\epsilon, \delta)$ DP for any $\epsilon$ and $\delta$ such that

$$ \epsilon \leq \min_{\alpha \in \mathbb{Z}, \alpha>1} \frac{\log \big(
\sum_{k=0}^{\alpha} \tbinom{\alpha}{k}(1-q)^{\alpha-k}q^k \mathcal{A}_{k}\big) + \log(1/\delta)}{1-\alpha}, $$

where $\mathcal{A}_{k} = \big(\mathbb{E}_{z \sim \mathsf{P}_0} \big[ (\mathsf{P}'_1(z)/{\mathsf{P}_0(z)})^k \big]\big)^{p}$.
Here, $\mathsf{P}_0 = \mathcal{N}(0,\sigma)*\mathcal{U}[-\tau/(2\eta),\tau/(2\eta)]$ and $\mathsf{P}'_1 = \mathcal{N}(c/\sqrt{p},\sigma)*\mathcal{U}[-\tau/(2\eta),\tau/(2\eta)]$, i.e., $\mathsf{P}_0$ shifted by $c/\sqrt{p}$.

See Appendix <ref>.

Corollary <ref> generalizes Theorem <ref> and recomputes the worst case divergence when we have both $l_2$ and $l_{\infty}$ norm sensitivity guarantees.
We prove that the worst case happens when the gradient $\nabla f(w_k,\tilde{x}, \tilde{y})$ of the differing element $(\tilde{x}, \tilde{y})$ is in a form whose Hamming weight is $p$ and each non-zero coordinate is equal to $\pm c/\sqrt{p}$.
Corollary <ref> also shows an efficient way to numerically compute the privacy loss by measuring the divergence between two one-dimensional distributions.

§.§ Privacy Amplification Examples

In Fig.  <ref>, we provide concrete examples of DP-SGD with ModelMix under various setups. We set $n=50,000$, $\delta=10^{-5}$, the clipping threshold $c=20$, and run DP-SGD under the Gaussian Mechanism <cit.>.
We measure the cumulative privacy loss in terms of $\epsilon$ under various sampling rates $q$, mixing thresholds $\tau$ and $l_{\infty}$ truncation parameter $p$. In all subfigures, we set the total number of iterations $T=5,000$, and the budget $\epsilon=200$ for regular DP-SGD.

In Fig. <ref> (a-c), we set the sampling rate $q = 0.02$ (a minibatch of size 1000 in expectation), and examine the privacy loss under $\tau = 0.075\eta, 0.15\eta, 0.3\eta$, where $\eta$ is the stepsize, respectively.
With ModelMix, especially with further help of Algorithm <ref> to ensure both $l_2$ and $l_{\infty}$ norm sensitivity, we achieve orders of magnitude improvement. For example, in Fig. <ref> (c) where $\tau = 0.3\eta$, compared to $(\epsilon=200, \delta=10^{-5})$ using regular DP-SGD, under the same setup, ModelMix, ModelMix with further $l_{\infty}$ truncation of $p=25$ and $p=100$ produce $\epsilon = 31.7, 5.4$ and $4.8$, respectively, with the same $\delta=10^{-5}$. The corresponding $\epsilon$ numbers produced in  Fig. <ref> (a) and (b) are (57.2, 17.9, 15.3) and (40.4, 9.0, 7.9), repectively.
The amplification factor of ModelMix matches our theoretical results, which is $\tilde{O}(1/\sqrt{\tau})$. On the other hand, the additional amplification from the $l_{\infty}$ norm truncation will reach some limit as $p$ increases. For large enough $p$, such combined amplification is empirically ${O}(1/{\tau})$. In Fig. <ref> (e-f), we change the sampling rate $q$ to be $0.04$ while keeping the remaining parameters the same as those in Fig. <ref> (a-c). In general, given larger privacy budget and mixing threshold $\tau$, ModelMix can produce stronger privacy amplification.

Privacy amplification from ModelMix on running DP-SGD with $n=50K$ samples under Gaussian Mechanism.

§.§ More Insights on Privacy Amplification

In the following, we provide more technical insights on how ModelMix improves the privacy. A quick intuition is that, compared to standard DP-SGD (Equation (<ref>)) where $w_k$ is only randomized by the noise $\Delta_{k}$, the sources of randomness in Algorithm <ref> are enriched, containing both the randomness in ModelMix and the noise perturbation.
More randomness often implies stronger privacy. However, though privacy relies on randomness, not all kinds of randomness can produce worst-case DP guarantees. ModelMix is an interesting example.
ModelMix must be applied with other noise mechanisms to produce meaningful DP guarantees, as explained below.

In ModelMix, the mixed state is within a bounded hull determined by earlier updates.
This means that the randomness is bounded and localized, which cannot produce reasonable worst-case DP guarantees, in contrast to Laplace/Gaussian noise.
To better illustrate this, consider the following example.
Suppose $\nabla F(w_k, \mathcal{D})=0$ and $\nabla F(w_k, \mathcal{D}')=1$ are two gradients evaluated by adjacent datasets $\mathcal{D}$ and $\mathcal{D}'$.
We randomize the gradients by introducing $\mathcal{U}[0,5]$, the uniform distribution between $[0,5]$.
The perturbed gradients now become equivalent to being uniformly sampled from $\mathcal{U}[0,5]$ and $\mathcal{U}[1,6]$ respectively.
When we compare the divergence between the two distributions $\mathcal{U}[0,5]$ and $\mathcal{U}[1,6]$, we cannot derive a worst-case $\epsilon$ bound ($\epsilon=\infty$ in this example) or an efficiently-composable $(\epsilon, \delta)$ DP guarantee [To achieve a meaningful privacy guarantee in DP-SGD for $T$ iterations, we need the failure probability $\delta = o(1/(nT))$ in a single iteration.]. To this end, additional Laplace or Gaussian noise is needed.

However, even with additional noise, we cannot improve the privacy loss of a single iteration significantly.
Continuing with the above example, in a single iteration of DP-SGD, if we add some Laplace noise $Lap(0,\lambda)$, which ensures $(\epsilon_0,0)$-DP, and further perturb the update by a uniform noise within $\mathcal{U}(0,5)$, the distributions of the two perturbed gradients $\nabla F(w_k, \mathcal{D})=0$ and $\nabla F(w_k, \mathcal{D}')=1$ now become $Lap(0,\lambda)*\mathcal{U}(0,5)$ and $Lap(0,\lambda)*\mathcal{U}(1,6)$, respectively.
Here, $*$ represents the convolution of two distributions. Still, the worst-case is not improved where we can only claim $\epsilon_0$-DP or still some not efficiently-composable $(\epsilon,\delta)$-DP for a single iteration.

So how does the localized randomness in ModelMix amplify the privacy? While ModelMix cannot significantly improve the worst-case utility-privacy tradeoff in a single iteration <cit.>, as shown in Theorems <ref> and <ref>, it smoothens the divergence between the output distributions from two arbitrary adjacent datasets.
Such an amplification is limited in a single iteration, but when we consider the composition of the privacy loss, the amplification accumulates to produce a strong improvement.
To be specific, for the DP-SGD update protocol (Equation (<ref>)), denoted by $\mathcal{M}$, we consider the  pointwise $\epsilon(w)$-loss at a particular output $w$ <cit.>,

$$ \epsilon_{\mathcal{D},\mathcal{D}',\mathsf{Aux}}(w) = \log \frac{\mathbb{P}(\mathcal{M}(\mathcal{D}, \mathsf{Aux})= w) }{\mathbb{P}(\mathcal{M}(\mathcal{D'},\mathsf{Aux})= w)},$$

where $\mathcal{D},\mathcal{D}'$ are two adjacent datasets and $\mathsf{Aux}$ represents the other auxiliary inputs used in $\mathcal{M}$.
When we take the supremum $\epsilon = \sup_{\mathcal{D},\mathcal{D}',\text{Aux}}\epsilon_{\mathcal{D},\mathcal{D}',\mathsf{Aux}}(w)$, it gives an equivalent $\epsilon$-DP guarantee. As explained before, incorporation of ModelMix  (Equation (<ref>)) cannot improve this (supremum) worst-case loss.
What we proved in Theorem <ref> is that, for any two datasets $\mathcal{D},\mathcal{D}'$ and $\mathsf{Aux}$, the expectation $\mathbb{E}[\epsilon(w)]$  and the variance of $\text{Var}[\epsilon(w)]$ for $w \sim \mathcal{M}(\mathcal{D},\text{Aux})$ will scale by $\tilde{O}(1/\tau)$ when ModelMix is further applied. Provided such an $\epsilon(w)$ loss of smaller mean and stronger concentration, we can derive a stronger high probability bound in a composite $(\epsilon, \delta)$ form when we measure the cumulative $\epsilon$ loss across the states $w_{[1:T]}$ from $T$ iterations.

§.§ Randomization beyond Noise

The analysis framework of ModelMix also sheds light on how to quantify the privacy amplification from a large class of “training-oriented” randomization commonly applied in deep learning.
There are many reasons to introduce randomness in optimization and learning other than privacy preservation.
For example, in stochastic gradient Langevin dynamics (SGLD) <cit.> for nonconvex optimization, it is common to utilize noisy gradient descent to escape saddle points <cit.>. Randomization can also strengthen the learning performance, e.g., random dropout <cit.> and data augmentation <cit.>.
In particular, data augmentation plays an important role in modern computer vision.
Generally speaking, data augmentation represents a large class of methods to improve robustness and reduce memorization (instead of generalization) by generating virtual samples through random cropping <cit.>, erasing <cit.> or mixing <cit.> the raw samples.

All of the above-mentioned randomnesses are localized, which cannot produce meaningful DP guarantees for the same reason as ModelMix.
For example, consider random erasing <cit.>, where a rectangle region of an image is randomly erased and replaced with random values.
When we process private images using the above mechanisms, the random transformation is multiplicative over the private input, meaning that the output is restricted to a bounded domain determined by the specific input processed.
Given two different images with at least one differentiating pixel, one can still distinguish them after random erasing as long as at least one differentiating pixel is not erased.
A similar argument holds for dropout, where a node in a neural network is ignored independently with some fixed rate, and the saddle point escaping algorithm <cit.>, where the gradient is perturbed by a bounded noise uniformly selected from a sphere. Our privacy analysis framework shows a way to analyze these randomizations in practical learning algorithms combined with DP-SGD to produce sharpened composition bounds.

§ FURTHER EXPERIMENTS

ModelMix is a generic technique, which can be implemented in almost all applications of DP-SGD without further assumptions.
In this section, we provide further experiments to measure its performance combined with other state-of-the-art advances in DP-SGD.
Below, when we report the improvement of ModelMix over existing works, we mostly follow the optimal hyper-parameters each work suggests.
We repeat each simulation five times and report the median of the results.
Details of the hyper-parameter selections we used can be found in Appendix <ref>. Based on the experiments conducted by previous works, to have a clear comparison, we will also test the proposed algorithms on the following three benchmark datasets, CIFAR10, SVHN and FMNIST[https://deepobs.readthedocs.io/en/stable/api/datasets/fmnist.html.].
CIFAR10 consists of 60,000 color images in 10 classes, where 50,000 are for training and 10,000 for test.
The Street View House Numbers (SVHN) dataset has 73,257 images of real world house digits for training and 26,032 for test.[We do not use the 600,000 auxiliary samples provided in SVHN in our experiments.]
Fashion MNIST (FMNIST) contains 70,000 greyscale images of fashion products from 10 categories with 60,000 for training and 10,000 for test.
In the following experiments, we will assume that all the training samples in the above-mentioned datasets are private.

§.§ Shallow Network

Since the magnitude of gradient perturbation is adversely dependent on the model size, instead of using the cutting-edge deep models, a large number of works are devoted to building small models to carefully balance model capacity and utility loss caused by DP-SGD <cit.>.
One of the best existing results is given by <cit.>.
In <cit.>, Tramer and Boneh showed that handcrafted features extracted from raw data can significantly strengthen the performance of training shallow models with DP-SGD.
They proposed to first privately use ScatterNet <cit.> to process the dataset and to then apply DP-SGD on the ScatterNet features.
With a DP budget of $(\epsilon=3,\delta=10^{-5})$,  <cit.> successfully trained a five-layer CNN with ScatterNet, which achieves $66.9\%$ and $87.2\%$ accuracy on CIFAR10 and FMNIST, respectively.
In the following, we consider further improving their results by using ModelMix.
With the same setup, we spend $\epsilon=0.695$ budget to privately estimate the statistics of datasets to apply ScatterNet.
In Table <ref>, we compare the utility-privacy tradeoff when we run DP-SGD with/without ModelMix on training the same CNN suggested by <cit.> on the ScatterNet features, where $\delta$ is always fixed to be $10^{-5}$ in all cases. Table <ref> shows that ModelMix can bring significant improvement even in simple model training with very low privacy budget.

Method$\backslash$Privacy     $\epsilon=0.2$     $0.4$     $0.6$     $0.8$     $1.0$    $ \infty$

<cit.>       32.8     45.9      56.1      61.2      64.1     74.6

DP-SGD +MM      58.5     63.3       65.0      66.9     68.3     74.6

Method$\backslash$Privacy     $\epsilon=0.2$     $0.4$     $0.6$     $0.8$     $1.0$     $\infty$

<cit.>      53.4     74.0      81.5      84.5      85.6     91.2

DP-SGD+MM      83.9     85.7       86.1      87.9     88.8     91.2

Test Accuracy of DP-SGD with/out ModelMix (MM) on training small CNN with ScatterNet Features.

§.§ Deep Models

In this subsection, we consider applying both ModelMix to help train large neural networks with DP-SGD. We take CIFAR10 and SVHN as examples. We will further apply Algorithm <ref> to enhance privacy amplification. In particular, according to our analysis on the clipping threshold and sampling noise, we select $c=20$ and the $l_{\infty}$-norm truncation parameter $p=100$ in all experiments.  In Table <ref>, we show the tradeoff between privacy and learning accuracy when training Resnet20 on CIFAR10 and SVHN with DP-SGD, respectively. Still, in all the cases, $\delta$ is fixed to be $10^{-5}$. With an $(\epsilon=8, \delta=10^{-5})$-DP budget, combined with ModelMix, we train Resnet20 which achieves $70.4\%$ and $90.1\%$ accuracy on CIFAR10 and SVHN, respectively. In comparison, regular DP-SGD can only produce $56.1\%$ and $84.9\%$ accuracy. In Table <ref>, we also include results for particular test accuracy the privacy budget required. In general, ModelMix roughly brings 20$\times$ improvement on $\epsilon$ to produce usable security parameters.

Method$\backslash$Privacy     $\epsilon=4$     $5$     $6$     $7$     $8$     $\infty$

Regular DP-SGD      47.3     51.1     53.2     54.7      56.1     90.7

DP-SGD+MM      59.8     64.6     67.5     69.3       70.4     90.7

Method$\backslash$Privacy     $\epsilon=2$     $4$     $5$     $6$     $8$     $\infty$

Regular DP-SGD      41.2      67.5      78.4     81.9     84.9     92.3

DP-SGD+MM      78.2     86.3     87.8     89.0     90.1       92.3

Test Accuracy of DP-SGD with/out ModelMix (MM) on training Resnet20.

Method$\backslash$Accuracy(%)     64     66     68     70     72

Regular DP-SGD      $\epsilon=95.8$     108.3     118.2     140.6      182.9

DP-SGD+MM      $\epsilon=4.7$     5.3     6.3       7.5     9.2

Method$\backslash$Accuracy(%)     $88$     $88.5$     $89$     $89.5$     $90.1$

Regular DP-SGD        $\epsilon = 93.6$      101.3     112.4     128.9     144.5

DP-SGD+MM      $\epsilon = 5.1 $     5.6     6.0     6.8       8

Privacy Loss of DP-SGD with/out ModelMix (MM) on training Resnet20.

§.§ Assistance with Public Data

With access to additional public data, many elegant ideas have been proposed to significantly improve the performance of DP-SGD. In this subsection, we present results where we use ModelMix to further improve two representative works <cit.> and <cit.>. <cit.> shows a private transfer learning method on CIFAR10, where a SimCLR model <cit.> is first pretrained on unlabeled ImageNet and a linear model is then trained on the features extracted from the
penultimate layer of the SimCLR model using DP-SGD. With a privacy budget $(\epsilon=2,\delta=10^{-5})$, one can achieve $92.7\%$ on CIFAR10 <cit.>. Similarly, we can apply ModelMix to improve this DP guarantee to $(\epsilon=0.64,\delta=10^{-5})$ with the same performance. On the other hand, if we assume the same privacy budget $(\epsilon=2,\delta=10^{-5})$ as <cit.>, we can instead achieve $93.6\%$, close to the non-private optimal performance $94.3\%$.

We also compare to <cit.>, which takes 2,000 ImageNet samples as public data and estimates a low-dimensional embedding of the private gradient when applying DP-SGD. Using the hyper-parameters suggested by <cit.>, we reproduce their experiments to privately train Resnet20 on CIFAR10, which achieves $73.2\%$ accuracy with an $(\epsilon=8,\delta=10^{-5})$ budget, and $79.1\%$ with $(\epsilon=111.2,\delta=10^{-5})$. With ModelMix, a sharpened tradeoff is produced, where with a budget $(\epsilon=2.9,\delta=10^{-5})$ and $(\epsilon=6.1,\delta=10^{-5})$, we achieve accuracy of $74.2\%$ and $79.1\%$, respectively. This is also close to the non-private optimal performance 82.3% of <cit.>, with the hyper-parameters suggested.

§ CONCLUSION AND PROSPECTS

In this paper, we present a formal study on the privacy amplification from the trajectory entropy and the influence of gradient clipping in DP-SGD.
We show fundamental improvements over DP-SGD, especially for deep learning, without assistance of other assumptions and additional data resources.
This is a first step to consider the amplification from the potential entropy underlying the intermediate computation in DP-SGD and there are still many interesting directions for further generalization of ModelMix.
For example, one may consider an adaptive selection of the envelope radius $\tau$ using techniques in <cit.> or approximate it with public data after a projection to a low-rank subspace <cit.>. Moreover, the privacy analysis presented can also be generalized to quantify amplification from a large class of practical randomness used in deep learning, as explained in Section <ref>; or to formalize the privacy guarantees of many heuristic privacy protections, such as Instahide <cit.> and Datamix <cit.>. Another issue we pointed out, which is of more interest to practitioners, is the important role of sampling noise in clipped DP-SGD and the connection to the clipping threshold. The theory and the empirical study shown here could be meaningful to instruct further architecture-level improvement, for example the application of batch norm for the noise variation reduction.  We leave this to future work.

§ PROOF OF THEOREM 3.1

For a $\beta$-smooth function $F(w)$, we have the following fact <cit.> that for any $w$ and $w'$,

\begin{equation}
\label{smooth_general}
-\frac{\beta}{2} \|w-w' \|^2   \leq F(w) - F(w') - \langle \nabla F(w'), w-w' \rangle \leq \frac{\beta}{2} \|w-w' \|^2.
\end{equation}

Equation (<ref>) will be constantly applied in the following proof.
It is noted that $\mathbb{E}[G_k]=nq\nabla F(w_k)$ and we use $\eta'$ to denote $\eta \cdot (nq)$ in the following.
Thus, $\mathbb{E}[\eta G_k] = \eta' \nabla F(w_k)$.

For each $k \in [0:T-1]$, conditional on $w_{k}$ and $w_{k-1}$, with the updating rule defined in Algorithm <ref>, we have

\begin{equation}
\label{GD-convergence}
\begin{aligned}
\mathbb{E}\big[ \|w_{k+1} - w^*\|^2\big]  \\
= &
\mathbb{E}\big[ \| {\alpha}_{k+1} w_{k} + (1-{\alpha}_{k+1}) w_{k-1}-w^* - \eta (G_{k} + \Delta_{k+1}) + e_{\tau_{k+1}} \|^2\big] \\
= &
\mathbb{E}\big[\| \alpha_{k+1} (w_{k}-w^*) + ({1}-{\alpha}_{k+1}) (w_{k-1}-w^*)\|^2 \\
- 2\eta' \langle \frac{w_{k}+w_{k-1}}{2}-w^*, \nabla F(w_{k})) \rangle \\
+ \frac{\eta'^2}{n^2q^2}\big(\mathbb{E}[\|G_{k}\|^2] + \mathbb{E}[\|\Delta_{k+1}\|^2]\big)  + \mathbb{E}[\|e_{\tau_{k+1}}\|^2].
\end{aligned}
\end{equation}

Here, $e_{\tau_{k+1}}$ represents the additional bias caused by the modification which ensures the per coordinate distance between $w_{k}$ and $w_{k-1}$ is at least $\tau_{k+1}$, whose mean is zero. Therefore, the variance $\mathbb{E}[\|e_{\tau_{k+1}}\|^2]$ of $e_{\tau_{k+1}}$ is bounded by $d\tau^2/12$.
Furthermore, we have that for each coordinate

\begin{equation*}
\begin{aligned}
\|{\alpha}_{k+1}(j) (w_{k}(j)-w^*(j)) + (1-{\alpha}_{k+1}(j)) (w_{k-1}(j)-w^*(j))\|^2\\
\leq
\alpha_{k+1}(j)\| w_{k}(j)-w^*(j)\|^2 + (1-{\alpha}_{k+1}(j)) \| w_{k-1}(j)-w^*(j)\|^2,
\end{aligned}
\end{equation*}

and thus,

\begin{equation*}
\begin{aligned}
\mathbb{E}[\| {\alpha}_{k+1} & (w_{k}-w^*) + ({1}-{\alpha}_{k+1}) (w_{k-1}-w^*)\|^2] \\
\leq
\frac{\|w_{k}-w^* \|^2 + \|w_{k-1}-w^* \|^2}{2}.
\end{aligned}
\end{equation*}

In addition, due to the convexity and the smoothness assumption,

\begin{equation}
\begin{aligned}
\langle  \frac{w_{k}+w_{k-1}}{2}-w^*,  \nabla F(w_{k}) \rangle \\
= & \langle w_{k}-w^* + \frac{w_{k-1}-w_{k}}{2} ,  \nabla F(w_{k}) - \nabla F(w^*) \rangle \\
\geq
F(w_{k}) - F(w^*) + \langle \frac{w_{k-1}-w_{k}}{2} ,  \nabla F(w_{k}) - \nabla F(w^*) \rangle\\
\geq & \frac{F(w_{k})+F(w_{k-1})}{2} - F(w^*) - \frac{\beta}{4}\|w_{k}-w_{k-1}\|^2.
\end{aligned}
\label{cross_bound}
\end{equation}

In (<ref>), we use (<ref>) and the following fact that $\langle w-w', \nabla F(w)-\nabla F(w') \rangle \geq F(w)-F(w')$ for any $w, w'$ and a convex function $F(\cdot)$.
With the Lipschitz assumption, we know $\|G_k\|^2 \leq n^2L^2$.

\begin{equation}
\begin{aligned}
& \mathbb{E}\big[\frac{F(w_{k})+F(w_{k-1})}{2} - F(w^*) \big] \\
\leq
\frac{1}{2\eta'} \big( \frac{\|w_{k}-w^*\|^2+\|w_{k-1}-w^*\|^2}{2} - \|w_{k+1}-w^*\|^2 + \frac{d\tau_{k+1}^2}{12} \big) \\
+ \frac{\eta'}{2n^2q^2} (n^2L^2 +  \mathbb{E}[\|\Delta_{k+1}\|^2]) + \frac{\beta}{4}\|w_{k}-w_{k-1}\|^2.
\end{aligned}
\label{GD-convergence-v3}
\end{equation}

Summing up both sides of Equation (<ref>) from $k=0,1,...,T-1$, we take expectation across $w_{[-1:T]}$ and have

\begin{equation*}
\label{GD-convergence-3}
\begin{aligned}
2 \mathbb{E} \big[ F(\frac{\sum_{k=1}^T w_{k-1}+w_{k-2}}{2T}) - F(w^*) \big] \\
\leq
{2\over T}\cdot  \sum_{k=1}^T  \mathbb{E}\big[\frac{F(w_{k-1})+F(w_{k-2})}{2} - F(w^*) \big] \\
\leq
\eta'^{-1} \big( \frac{\|w_{-1}-w^*\|^2+2\|w_{0}-w^*\|^2}{2T} +  \frac{\sum_{k=0}^{T-1}d\tau^2_{k+1}}{12T} \big) + \frac{\eta'L^2}{q^2} \\
& + \frac{\eta'}{n^2q^2} ( \mathbb{E}[\sum_{k=0}^{T-1}\|\Delta_{k+1}\|^2/T]) + \frac{\beta\mathbb{E}[\sum_{k=0}^{T-1}\|w_{k}-w_{k-1}\|^2]}{2T}.\\
\end{aligned}
\end{equation*}

It is noted that the updating rule can be written as,

\begin{equation}
\begin{aligned}
w_{k+1} - w_{k}
& = -(1-\alpha_{k+1}) (w_{k}-w_{k-1}) - \eta(G_k + \Delta_{k+1}) \\
& = (1-\alpha_{k+1})(1-\alpha_{k})(w_{k-1}-w_{k-2})+ \\
& ~~~ (1-\alpha_{k+1})\eta(G_{k-1} + \Delta_{k}) - \eta(G_k + \Delta_{k+1})
\end{aligned}
\end{equation}

By the recursion and the independce of different $\alpha_{k}$, we have a closed-form upper bound of $\| w_{k+1} - w_{k} \|^2$ conditional on the initialization $\|w_{0}-w_{-1}\|^2$, where

\begin{equation}
\begin{aligned}
\label{difference recursion}
& ~~~~\mathbb{E}[\|w_{k+1} - w_{k}\|^2] \\
&\leq(\frac{1}{3})^{k+1}\|w_{0}-w_{-1}\|^2  + \frac{11}{2}\eta^2(n^2L^2+\mathbb{E}[\|\Delta\|^2]) \\
& ~~~~+ 2(k+1)(1/2)^{k+1}\eta\|w_{0}-w_{-1}\|(nL+\mathbb{E}[\|\Delta\|]).
\end{aligned}
\end{equation}

Here, in (<ref>), we use the following fact that

\[
\mathbb{E}[\|\sum_{i=1}^{k} \prod_{j=1}^{i-1} \alpha_{j}\|^2] = \sum_{i=1}^k (\frac{1}{3})^i + 2\sum_{i=1}^k \sum_{j=1}^{i} (\frac{1}{2})^j(\frac{1}{3})^{i-j} \leq \frac{3}{2} + 4 = \frac{11}{2},
\]

Now putting all the above together, we have that $\mathbb{E} \big[  F(\frac{\sum_{k=1}^T w_{k-1}+w_{k-2}}{2T}) - F(w^*) \big]$ is upper bounded by

\[
\small{
\begin{aligned}
&\frac{\|w_{-1}-w^*\|^2+2\|w_{0}-w^*\|^2 + \sum_{k=1}^Td\tau^2_{i}/12}{2\gamma \sqrt{T}} \\
& + \frac{\gamma(L^2/q^2+\mathbb{E}[\|\Delta\|^2]/(n^2q^2))}{2\sqrt{T}} + \frac{\beta\big( \frac{11}{2}\eta^2(n^2L^2+\mathbb{E}[\|\Delta\|^2]) \big)}{T}\\
& + \frac{\beta \big(\frac{3}{2}\|w_0 - w_{-1}\|^2 + 4\eta \|w_{0}-w_{-1}\|(nL+\mathbb{E}[\|\Delta\|])}{4T} \\
& = \frac{\|w_{-1}-w^*\|^2+2\|w_{0}-w^*\|^2 + \sum_{k=1}^T d \tau^2_{i}/12}{2\gamma \sqrt{T}}\\
& ~~~~+  \frac{3\beta\|w_0 - w_{-1}\|^2+ 11\gamma^2(\frac{L^2}{q^2}+\mathbb{E}[\frac{\|\Delta\|^2}{n^2q^2}])}{8T} \\
& ~~~~+ \frac{\beta\gamma \|w_{0}-w_{-1}\|(\frac{L}{q}+\mathbb{E}[\frac{\|\Delta\|}{nq}])}{T^{3/2}}+ \frac{\gamma(L^2+\mathbb{E}[\|\Delta\|^2])}{2\sqrt{T}}.
\end{aligned}}
\]

§ PROOF OF THEOREM 3.2

We first prove the following theorem, which analyzes the standard clipped DP-SGD without ModelMix.

Suppose the objective loss function $F(w)$ is $\beta$-smooth and satisfies Assumption <ref>, then there exists some constant $\psi > 0$ such that when the clipping threshold $c$ satisfies

\begin{equation}
\label{c_req}
c \geq \max \{4\kappa\log(10), -\psi\kappa\log(\kappa)\log(\frac{\sqrt{d \log(1/\delta)}}{n\epsilon })\},
\end{equation}

then the convergence rate of per-sample clipped DP-SGD with an $(\epsilon, \delta)$-DP guarantee satisfies

\begin{equation}
\begin{aligned}
\label{batch_main}
\mathbb{E} \Big[ &\frac{\sum_{k=0}^{T-1}  \min\big\{ 9/20 \cdot  \|\nabla F(w_k)\|^2, c/20 \cdot \| \nabla F(w_k)\| \big \}}{T}\Big]  \\
&\leq (\frac{v}{2} +\frac{5}{2})\cdot \frac{c\sqrt{\mathcal{R}_{F}\beta d \log(1/\delta)}}{n\epsilon},
\end{aligned}
\end{equation}

where $\mathcal{R}_{F} = \sup_{w}F(w)- \inf_{w} F(w)$ and $v$ is some constant determined by noise Mechanism.

At the $k$-th iteration, based on the smoothness assumption, we have

\begin{equation}\label{equ:smooth}
F(w_{k+1}) \leq F(w_{k}) + \langle \nabla F(w_{k}), w_{k+1}-w_{k}\rangle + {\beta\over 2} \cdot \|w_{k+1}-w_{k} \|^2.
\end{equation}

We consider the following equivalent updating rule, where the step size $\eta$ is scaled by $1/(nq)$ in  (<ref>),

\[
w_{k+1} - w_{k} = \eta \cdot \Big( \frac{1}{nq} \cdot (\sum_{j=1}^{B_k} \mathsf{CP}(g^j_k, c) + \Delta_{k+1}) \Big).
\]

Here $g^j_k = \nabla f(w_k, x_{(j)}, y_{(j)})$, where $(x_{(j)}, y_{(j)})$ is the $j$-th sample selected in the minibatch $S_k$. For simplicity, we use ${g}_{k}$ to denote the random variable

\[
g_k =  \nabla f(w_k, x, y),
\]

where $(x, y)$ is randomly selected from the dataset $\mathcal{D}$. Thus, we have the following observation that

$$ \mathbb{E}[\mathsf{CP}(g_k, c)] = \mathbb{E}[\frac{1}{nq} \cdot \sum_{j=1}^{B_k} \mathsf{CP}(g^j_k, c)],$$

due to the i.i.d. sampling of rate $q$.  Thus, with expectation we have

\[
\mathbb{E}[w_{k + 1} - w_k]
= \eta \Big( \frac{1}{B_k} \sum_{j=1}^{B_k} \mathbb{E}[\mathsf{CP}(g^j_k, c)] + \mathbb{E}[\Delta_{k+1}] \Big)
= \eta \mathbb{E}[\mathsf{CP}(g_k, c)],
\]

\[
\mathbb{E}[\|w_{k+1}-w_{k} \|^2]
\leq
\eta^2 \cdot (c^2/q^2 + \mathbb{E}[\|\Delta_{k+1}\|^2]/(n^2q^2)).
\]

Taking this into Equation (<ref>), we have

\begin{equation}
\label{smooth_basic}
\begin{aligned}
\mathbb{E}[F(w_{k+1})]
& \leq F(w_k) - \mathbb{E}[\langle \nabla F(w_{k}), w_{k + 1} - w_k \rangle] \\
&   \quad + {\beta\over 2} \cdot \mathbb{E}[\|w_{k+1}-w_{k} \|^2] \\
& \leq F(w_k) - \eta \cdot\mathbb{E}[\langle \nabla F(w_{k}), \mathsf{CP}(g_{k}, c) \rangle] \\
&   \quad + {\beta \eta^2\over 2} \cdot (c^2/q^2 + \mathbb{E}[\|\Delta_{k+1}\|^2/(n^2q^2)]).
\end{aligned}
\end{equation}

Let  $\gamma_k = \min \{1, {c\over \|g_k\|}\}$, i.e., $\mathsf{CP}(g_{k}, c)= \gamma_k g_k$.
From previous works <cit.> on the standard Gaussian DP mechanism, we know that the variance of the noise $\Delta_k$ is bounded by $O(q^2dT\log(1/\delta) / \epsilon^2)$.
Therefore, there must exist some constant $v$ such that we can rewrite Equation (<ref>) as

\begin{equation}
\label{smooth_basic_1}
\begin{aligned}
\mathbb{E}[F(w_{k+1})]
& \le F(w_k) - \eta \cdot \mathbb{E}[\langle \nabla F(w_{k}), \gamma_k g_{k} \rangle] \\
& \quad + {\beta \eta^2\over 2} \cdot (\mathbb{E}[\|\gamma_k  g_k \|^2] + \mathbb{E}[\|\Delta_{k+1}]\|^2) \\
& \leq  F(w_k) -  \eta\cdot \mathbb{E}[\langle \nabla F(w_{k}), \gamma_k g_{k} \rangle] \\
& \quad + {\beta \eta^2\over 2}\cdot (c^2/q^2 + v\cdot \frac{dT\log(1/\delta)}{n^2\epsilon^2}).
\end{aligned}
\end{equation}

In the following Lemma, we lower bound $\mathbb{E}[\langle \nabla F(w_{k}), \gamma_k  g_{k} \rangle]$. We use $\xi_k = g_k - \nabla F(x_k)$ to denote the stochastic gradient noise.

For any $p > 2$ and $c_0 > 0$ such that $c= p\cdot c_0$, it holds that

$$\mathbb{E}[\langle \nabla F(w_{k}), \gamma_k  g_{k} \rangle] \geq \min \{a_1\|\nabla F(x_k) \|^2 - a_2,  a_3\|\nabla F(x_k)\| \},$$

where $a_1$, $a_2$ and $a_3$ are constants define as below

\[
a_1 = \frac{(1-e^{-{c_0\over \kappa}})}{2},
a_2 = \Big(\frac{(c_0+\kappa)\cdot e^{-{c_0\over\kappa}}}{1-e^{-{c_0\over\kappa}}}\Big)^2,
\]

\[
a_3 = c\Big( (1-e^{-{c_0\over\kappa}})(1-\frac{1}{p-1}-\frac{1}{p-2})- e^{-{c_0\over\kappa}} \Big).
\]

The proof contains two parts.

(1). First, we consider the case when the gradient is small where $\|\nabla F(x_k)\| \leq (p-1)c_0$.

\begin{equation}
\label{small_gradient}
\begin{aligned}
\mathbb{E}[\langle \nabla F(w_{k}), \gamma_k  g_{k} \rangle] \\
& = \mathbb{E}[\langle \nabla F(w_{k}), \gamma_k (\nabla F(w_{k})+\xi_k ) \rangle] \\
& = \|\nabla F(w_{k})\|^2 \cdot \mathbb{E}[\gamma_k] + \mathbb{E}[\langle \nabla F(w_k) ,\xi_k\rangle \cdot (1-\gamma_k)] \\
& \geq \|\nabla F(w_{k})\|^2 \Pr(\|g_k\| \leq c) - \| \nabla F(w_k)\|\mathbb{E}[\|\xi_k\|\cdot \bm{1}_{\|g_k\|>c}] \\
& \geq \|\nabla F(w_{k})\|^2 \Pr(\|\xi_k\| \leq c_0) - \| \nabla F(w_k)\|\mathbb{E}[\|\xi_k\|\cdot \bm{1}_{\|g_k\|>c}]
\end{aligned}
\end{equation}

Here, the second equality is because $\mathbb{E}[\langle \nabla F(w_k), \xi_k]=0$. In the third inequality, we use the following facts.
It is not hard to see that

\[
\|\nabla F(w_{k})\|^2 \cdot \mathbb{E}[\gamma_k] \geq \|\nabla F(w_{k})\|^2 \cdot \Pr(\|g_k\|\leq c).
\]

Since $\|F(w_k)\| \leq (p-1)c_0$ is assumed, we have the probability of clipping

\[
\Pr(\|g_k\|\leq c) \geq \Pr(\|\xi_k\|\leq c_0).
\]

As for the second term, when $\|g_k\|\leq c$, i.e., no clipping is performed.
In this case, $\gamma_k = 1$ and we have

\[
\mathbb{E}[\langle \nabla F(w_k) ,\xi_k\rangle \cdot (1-\gamma_k) \cdot \bm{1}_{\|g_k\|\leq c}] = 0.
\]

Thus, we only need to consider $\mathbb{E}[\langle \nabla F(w_k) ,\xi_k\rangle \cdot (1-\gamma_k) \cdot \bm{1}_{\|g_k\|\leq c}]$ when clipping is performed.
By Cauchy-Schwartz inequality and $1-\gamma_k \leq 1$, we have the bound in Equation (<ref>). Now, we turn to upper bound $\mathbb{E}[\|\xi_k\|\cdot \bm{1}_{\|g_k\|>c}] $. By the Integrated Tail Probability Expectation Formula, we can derive the following variant,

\begin{equation*}
\begin{aligned}
\mathbb{E}[\|\xi_k\|\bm{1}_{\|g_k\|>c}]
& \leq \mathbb{E}[\|\xi_k\|\bm{1}_{\|\xi_k\|>c_0}] \\
& = c_0(1-\Pr(\|\xi_k\|<c_0)) + \int_{c_0}^{+\infty} \Pr(\|\xi_k\|>t) dt  \\
& \leq c_0 e^{-c_0/\kappa} + \int_{c_0}^{+\infty} e^{-t/\kappa} dt  = (c_0 + \kappa) e^{-c_0/\kappa}.
\end{aligned}
\end{equation*}

Putting things together, we have

\begin{equation*}
\begin{aligned}
\mathbb{E}[\langle \nabla F(w_{k}), \gamma_k  g_{k} \rangle] \\
& \geq (1-e^{-c_0/\kappa})\|\nabla F(w_k)\|^2 -(c_0 + \kappa) e^{-c_0/\kappa} \|\nabla F(x_k)\| \\
& \geq \frac{(1-e^{-c_0/\kappa})}{2}\|\nabla F(w_k)\|^2  - (\frac{(c_0+\kappa)e^{-c_0/\kappa}}{1-e^{-c_0/\kappa}})^2.
\end{aligned}
\end{equation*}

Here, the second inequality we use the following fact: $ax^2-bx \geq (a/2) x^2 -(b/a)^2,$ for any positive constants $a$.

(2). Second, we consider the case when the gradient is larger  where $\|\nabla F(x_k) \|>(p-1)c_0$.

\begin{equation}
\label{large_gradient}
\small{
\begin{aligned}
\mathbb{E}[\langle \nabla F(w_{k}), \gamma_k  g_{k} \rangle] \\
= &
\|\nabla F(w_{k})\|^2 \cdot \mathbb{E}[\gamma_k] + \mathbb{E}[\gamma_k\langle \nabla F(w_k) ,\xi_k\rangle] \\
\geq &
\|\nabla F(w_{k})\|^2 \mathbb{E}[\gamma_k  \cdot \bm{1}_{\|\xi_k\|\leq c_0}] - \| \nabla F(w_k)\|\mathbb{E}[\gamma_k \|\xi_k\|\cdot \bm{1}_{\|\xi_k\|\leq c_0}] \\
+ \mathbb{E}[ \langle \nabla F(w_k),\gamma_k g_k\rangle\cdot \bm{1}_{\|\xi_k\|>c_0}] \\
\geq &
c \|\nabla F(w_k)\|\Pr(\|\xi_k\| \leq c_0) \cdot \big ( \frac{\| \nabla F(w_k)\| }{\| \nabla F(w_k)\|+c_0} - \frac{c_0}{\|\nabla F(w_k)\|-c_0}  \big ) \\
+ \mathbb{E}[ \langle \nabla F(w_k),\gamma_k g_k\rangle\cdot \bm{1}_{\|\xi_k\|>c_0}]. \\
\end{aligned}
\end{equation}

Here, the second inequality is because the term

\[
\big ( \frac{\| \nabla F(w_k)\|}{\|\nabla F(w_k)\| + \|\xi_k\|} - \frac{\|\xi_k \|}{\|\nabla F(w_k)\|-\|\xi_k\|}  \big )
\]

reaches the minimal when  $\|x_{i_k}\|=c_0$ for $\|x_{i_k}\| \in [0,c_0]$.
On one hand,

$$  \frac{\| \nabla F(w_k)\| }{\| \nabla F(w_k)\|+c_0} - \frac{c_0}{\|\nabla F(w_k)\|-c_0}  \geq  \frac{(p-1)}{p} - \frac{1}{p-2}.$$

Moreover, we have

$$ \mathbb{E}[ \langle \nabla F(w_k),\gamma_k g_k\rangle\cdot \bm{1}_{\|\xi_k\|>c_0}] \geq -c\| \nabla F(w_k)\| \Pr(\|\xi_k\|>c_0),$$

since $\|\gamma_k g_k\|<c$.
Thus, with the assumption $\Pr(\|\xi_k\|> c_0) \leq e^{-c_0/\kappa}$, we have

\begin{equation*}
\begin{aligned}
&\mathbb{E}[\langle \nabla F(w_{k}), \gamma_k  g_{k} \rangle]  \\
\geq &
c\|\nabla F(w_{k})\|\big( (1-e^{-c_0/\kappa})(1-\frac{1}{p-1}-\frac{1}{p-2})- e^{-c_0/\kappa} \big).
\end{aligned}
\end{equation*}

Now, we pick $p=4$, and we have the following corollary, i.e.,

\begin{equation*}
\mathbb{E}[\langle \nabla F(w_{k}), \gamma_k  g_{k} \rangle] + a_2 \geq \min \{a_1\|\nabla F(x_k) \|^2,  a_3\|\nabla F(x_k)\| \},
\end{equation*}

\[
a_1 = \frac{(1-e^{-c/(4\kappa)})}{2},
a_2 = (\frac{(c/4+\kappa)e^{-c/(4\kappa)}}{1-e^{-c/(4\kappa)}})^2
\]

\[
a_3 = c\big(1/6 -7/6 \cdot e^{-c/(4\kappa)}) \big).
\]

Now, we return to Equation (<ref>) and sum up across $k=0,1,...,T-1$, and we have

\begin{equation}
\label{dp_tradeoff}
\begin{aligned}
\mathbb{E}[ \frac{\sum_{k=0}^{T-1} \min\big\{ a_1 \|\nabla F(w_k)\|^2, a_3\| \nabla F(w_k)\| \big \}}{T}] \\
& \leq \frac{F(w_0)-F(w_T)}{T\eta} + \frac{\beta \eta}{2} \big( \frac{c^2}{q^2} + v\frac{dT\log(1/\delta)}{n^2 \epsilon^2} \big) + a_2 \\
& = (\frac{v}{2} +\frac{3}{2})\frac{c\sqrt{\mathcal{R}_{F}\beta d \log(1/\delta)}}{n\epsilon q^2} + a_2.
\end{aligned}
\end{equation}

Here, we set $\eta = \frac{\sqrt{\mathcal{R}_{F}d \log(1/\delta)}}{n \epsilon (c/q) \sqrt{\beta}}$ and $T = \frac{(n\epsilon)^2}{d\log(1/\delta)}$. Finally, we want to pick $c$ to ensure $a_2$ is at most in the same order of $O(\frac{\sqrt{d\log(1/\delta)}}{n\epsilon})$, while $a_1$ and $a_3$ are positive.
To this end, we select $c$ such that

$$ c \geq  4\kappa\log(10),
\frac{c}{2\kappa} - \log\frac{(c/4+\kappa)^2}{c} \geq 0.4 - \log(\frac{\sqrt{\mathcal{R}_{F}\beta d \log(1/\delta)}}{n\epsilon}), $$

then the right hand of (<ref>) is further bounded by $(\frac{v}{2} +\frac{5}{2})\frac{c\sqrt{\mathcal{R}_{F}\beta d \log(1/\delta)}}{n\epsilon}$ while $a_1 = 9/20$ and $a_3 = c/20.$

We still consider the case where the step size $\eta$ is scaled by a factor $1/(nq)$ for simplicity.
For each $k \in [1:T]$, from the updating rule in Algorithm <ref>, we have that conditional on $w_{k}$ and $w_{k-1}$,

\begin{equation}
\label{GD-convergence-1}
\footnotesize{
\begin{aligned}
& \mathbb{E}\big[ F(w_{k+1}) \big]  \\
& \leq \mathbb{E}\big [F(w_{k}) + \langle \nabla F(w^{k}), w^{k+1}-w^{k}\rangle + \frac{\beta}{2}\|w^{k+1} - w^{k} \|^2 \big] \\
& = F(w_{k}) + \langle \nabla F(w_k), \frac{w_{k-1}-w_{k}}{2} - \eta \mathbb{E}[\mathsf{CP}(g_k, c)] \rangle \\
& + \frac{\beta}{2} \mathbb{E}\big[\|(1-\alpha_{k+1})(w_{k-1}-w_{k}) - \eta/(np)(G_k  + \Delta_{k+1}) + e_{\tau_{k+1}} \|^2 \big] \\
& \leq  F(w_{k}) + \frac{\langle \nabla F(w_{k}),w_{k-1}-w_{k}  \rangle }{2} -\eta \langle \nabla F(w_k), \mathbb{E}[\mathsf{CP}(g_k, c)]  \rangle \\
& \quad + \beta\big( \frac{\|w_{k}-w_{k-1} \|^2}{3}
+ \eta^2(c^2/q^2 + \|\Delta_{k+1}\|^2/(n^2q^2))+\|e_{\tau_{k+1}}\|^2 ] \big) .
\end{aligned}
\end{equation}

Here, we use the fact that $\mathbb{E}(1-\alpha)^2 = 1/3$ for a random $\alpha \in (0,1)$ and the AM-GM inequality. Still, $e_{\tau_{k+1}}$ represents the additional bias caused by the modification which ensures the per coordinate distance between $w_{k}$ and $w_{k-1}$ is at least $\tau_{k+1}$. We have  that the variance $\mathbb{E}[\|e_{\tau_{k+1}}\|^2]$ of $e_{\tau}$ is bounded by $d\tau_{k+1}^2/12$.

Now we apply (<ref>) again on the term $\frac{\langle \nabla F(w^{k}),w_{k}-w_{k-1}  \rangle }{2}$, we have that

\begin{equation}
\frac{\langle \nabla F(w^{k}),w_{k-1}-w_{k}  \rangle }{2} \leq \frac{1}{2} \big ( F(w_{k-1})- F(w_{k}) + \frac{\beta}{2}\|w_{k} - w_{k-1}\|^2\big).
\label{reverse_smooth}
\end{equation}

Now, substitute (<ref>) back to (<ref>), we have

\begin{equation}
\begin{aligned}
& \mathbb{E}[F(w_{k+1})]\leq \frac{F(w_k)+F(w_{k-1})}{2}  -\eta \langle \nabla F(w_k), \mathbb{E}[\mathsf{CP}(g_k, c)]  \rangle \\
& \quad + \beta\big(\frac{7 \| w_k-w_{k-1}\|^2 + \tau^2_{k+1}}{12} + \frac{\eta^2}{q^2} \cdot  (c^2+ \mathbb{E}[\|\Delta_{k+1}\|^2]/n^2) \big).
\end{aligned}
\label{GD-convergence-2}
\end{equation}

In the above inequality, we use the fact that $\|\mathsf{CP}(g^j_k,c)\| \leq c$ due to the clipping. We sum up both sides of (<ref>) and have

\begin{equation}
\begin{aligned}
& \mathbb{E}\big[\frac{\eta \sum_{k=0}^{T-1} \langle \nabla F(w_k), \mathsf{CP}(g_k,c) \rangle}{T}\big]  \leq \frac{3\mathcal{R}_{F}}{2T}+\beta \eta^2(\frac{c^2}{q^2}+\frac{\mathbb{E}[\|\Delta\|^2}{(nq)^2})\\
& + \frac{\beta\sum_{k=0}^{T-1}\big(
d\tau^2_{k+1} + 7\|w_{k}-w_{k-1} \|^2 \big)}{12T}
\end{aligned}
\label{MM_CSGD-1}
\end{equation}

where $\mathcal{R}_{F} = \sup_{w}F(w)- \inf_{w} F(w)$.
With a similar reasoning as  (<ref>) in Appendix <ref>, we have that given $w_{0}$ and $w_{-1}$,

\begin{equation}
\begin{aligned}
\sum_{k=0}^{T-1} \| w_k-w_{k-1}\|^2 & \leq \frac{3}{2} \|w_0-w_{-1}\|^2 + 4\eta c\|w_0-w_{-1} \| \\
& + \frac{11T}{2}\eta^2(c^2/q^2+\mathbb{E}[\|\Delta\|^2]/(nq)^2).
\end{aligned}
\label{additional_term_smooth}
\end{equation}

The rest of the analysis for the term  $\mathbb{E}\big[\frac{\eta \sum_{k=0}^{T-1} \langle \nabla F(w_k), \mathsf{CP}(g_k,c) \rangle}{T}\big]$ is the same as the proof of Theorem <ref> where virtually we handle a function whose smooth parameter becomes $\frac{101}{12}\beta$. Therefore, we still select
$\eta = \frac{\sqrt{\mathcal{R}_{F}d \log(1/\delta)}}{n \epsilon (c/q) \sqrt{\frac{101}{12}\beta}}$ and $T = \frac{(n\epsilon)^2}{d\log(1/\delta)}$, and substitute them into (<ref>), and then the right hand of (<ref>) in the
ModelMix case becomes

\begin{equation}
\small{
\begin{aligned}
(\frac{v}{2} +\frac{5}{2})&\frac{c\sqrt{\mathcal{R}_{F}\frac{101}{12}\beta d \log(1/\delta)}}{n\epsilon} + \frac{28c\beta{d\log(1/\delta)}\|w_0-w_{-1} \| }{12 q(n\epsilon)^2} \\
&  +  \frac{c d\log(1/\delta)\sqrt{\frac{101}{12}}\beta^{3/2}}{qn\epsilon\sqrt{\mathcal{R}_{F}}} \big(\frac{\sum_{k=1}^T d\tau^2_{k}}{12} + \frac{21\|w_0-w_{-1}\|^2}{24}\big) .
\end{aligned}
\end{equation}

§ PROOF OF THEOREM <REF> AND COROLLARY <REF>

For simplicity, we normalize $\bar{\tau} = \tau/\eta$. We consider the following equivalent aggregation model $\mathcal{M}$. Let $\mathcal{D}=(a_1,a_2,...,a_n)$ and $\mathcal{D}'=(a_1,a_2,...,a_{n-1})$ be two adjacent datasets. $\mathcal{M}(\mathcal{D})$ ($\mathcal{M}(\mathcal{D}')$) implements as follows. First, we apply i.i.d. sampling of rate $q$ to generate a subset $S$ from $\mathcal{D}$ $(\mathcal{D}')$ and output its sum $\Sigma(S) = \sum_{a_i \in S} a_i$ perturbed by a sum of Gaussian and independent ModelMix, i.e., $\mathcal{U}[-\frac{\bar{\tau}}{2}, \frac{\bar{\tau}}{2}]^d \times \mathcal{N}(\bm{0}, \sigma^2 \cdot \bm{I}_{d}),$ where we denote this kind of mixture distribution as $\mathcal{MG}^d_{\bar{\tau},\sigma}(\bm{0})$.  In the above model, $a_i = \mathsf{CP}(\nabla f(w,x_i,y_i), c)$ represents the individual gradient of each datapoint clipped by c for any $i$, where we model the distribution of DP-SGD with ModelMix procedure equivalently as

\begin{equation}
\mathcal{U}^{d}[-\frac{\tau}{2},\frac{\tau}{2}] + \eta(\sum_{ a_i \in S} a_i + \Delta).
\label{agg_model}
\end{equation}

Now, we can rewrite the distributions of $\mathcal{M}(\mathcal{D})$ and $\mathcal{M}(\mathcal{D}')$ as follows,

$$ \mathcal{M}(\mathcal{D}') = \sum_{S} p_S \mathcal{MG}^d_{\bar{\tau},\sigma}(\Sigma(S)),$$

for any $S \subset \mathcal{D} \cap \mathcal{D}'$ and $P_S$ represents the probability that $S$ gets sampled from $\mathcal{D}'$ with i.i.d. sampling of rate $q$. Correspondingly,

$$  \mathcal{M}(\mathcal{D}) = \sum_{S} p_S \big( (1-q)\mathcal{MG}^d_{\bar{\tau},\sigma}(\Sigma(S)) + q \mathcal{MG}^d_{\bar{\tau},\sigma}(\Sigma(S)+ a_n) \big). $$
Therefore, when we measure the $\alpha$-Rényi divergence, we have

\begin{equation}
\label{Rényi-shift}
\begin{aligned}
& \mathsf{D}_{\alpha}(\mathcal{M}(\mathcal{D})\|\mathcal{M}(\mathcal{D}'))\\
&\leq \max_{S} \mathsf{D}_{\alpha} \big(\mathcal{MG}^d_{\bar{\tau},\sigma}(\Sigma(S))\|  (1-q)\mathcal{MG}^d_{\bar{\tau},\sigma}(\Sigma(S)) + q \mathcal{MG}^d_{\bar{\tau},\sigma}(\Sigma(S, a_n) ) \big) \\
& \leq \max_{v,\|v\|\leq c} \mathsf{D}_{\alpha} \big(\mathcal{MG}^d_{\bar{\tau},\sigma}(\bm{0})\|  (1-q)\mathcal{MG}^d_{\bar{\tau},\sigma}(\bm{0}) + q \mathcal{MG}^d_{\bar{\tau},\sigma}(v) \big).
\end{aligned}
\end{equation}

In the first and the second inequality of (<ref>), we use the quasi-convexity and translation invariance properties of Rényi divergence, respectively. Thus, in the second inequality we subtract $\Sigma(S)$ on both side and we know $\|v = \Sigma(S \cup a_n) - \Sigma(S) = a_n\|$ is upper bounded by $c$.

To proceed, we first use the result (Theorem 5) in <cit.>, which suggests that $\mathsf{D}_{\alpha}(\mathcal{M}(\mathcal{D})\|\mathcal{M}(\mathcal{D}')) \geq \mathsf{D}_{\alpha}(\mathcal{M}(\mathcal{D}')\|\mathcal{M}(\mathcal{D}))$ and it suffices to consider $\mathsf{D}_{\alpha}(\mathcal{M}(\mathcal{D})\|\mathcal{M}(\mathcal{D}'))$ to derive the RDP bound in the following. It is noted that

\[
\small{
\begin{aligned}
& (1-\alpha)\mathsf{D}_{\alpha}(\mathcal{M}(\mathcal{D})\|\mathcal{M}(\mathcal{D}')) \\
& \leq \log \max_{v} \mathbb{E}_{o \sim \mathcal{MG}^d_{\tau,\sigma}(\bm{0})} \big( (1-q) + q\frac{\mathcal{MG}^d_{\bar{\tau},\sigma}(v)(o)}{\mathcal{MG}^d_{\bar{\tau},\sigma}(\bm{0})(o)} \big)^{\alpha}\\
& = \log \max_{v} \big\{\sum_{k=0}^{\alpha} \tbinom{\alpha}{k}(1-q)^{\alpha-k}q^k \mathbb{E}_{o \sim \mathcal{MG}^d_{\bar{\tau},\sigma}(\bm{0})}[\big(\frac{\mathcal{MG}^d_{\bar{\tau},\sigma}(v)(o)}{\mathcal{MG}^d_{\bar{\tau},\sigma}(\bm{0})(o)}\big)^{k}] \big\}.
\end{aligned}
\]

Thus, in the following, it is equivalent to consider $\mathcal{A}_k =\mathbb{E}_{o \sim \mathcal{MG}^d_{\bar{\tau},\sigma}(\bm{0})}[\big(\frac{\mathcal{MG}^d_{\tau,\sigma}(v)(o)}{\mathcal{MG}^d_{\bar{\tau},\sigma}(\bm{0})(o)}\big)^{k}],$ which has a semi-closed form

\begin{equation}\label{moment_mm}
\small{
\begin{aligned}
\mathcal{A}_k & = \max_{v} \sum_{j=1}^d \log \int  \frac{\big(\mathbb{P}(\mathcal{MG}^1_{\bar{\tau},\sigma}(v_j)={o}_j)\big)^{k}}{\big(\mathbb{P}(\mathcal{MG}^1_{\bar{\tau},\sigma}(0) ={o}_j)\big)^{k-1}} d o_j = \max_{v} \sum_{j=1}^d h(v_j).
\end{aligned}
\end{equation}

Here, $v_j$ and $o_j$ represents the $j$-th coordinated of $v$ and $o$, and we use the fact that the distribution of each coordinate of $\mathcal{MG}^d_{\bar{\tau},\sigma}(v)$ or $\mathcal{MG}^d_{\bar{\tau},\sigma}(\bm{0})$ is independent and thus its density function is a product of $\mathcal{MG}^1_{\bar{\tau},\sigma}(v_j)$ or $\mathcal{MG}^1_{\bar{\tau},\sigma}(0)$. To be specific, we have that $h(v_j)$ can be expressed as

\[
h(v_j) = \log \int_{o}\frac{\big(\int_{-\frac{\bar{\tau}}{2}}^{\frac{\bar{\tau}}{2}}\frac{e^{-(o-\alpha-v_j)^2/(2\sigma^2)}}{\bar{\tau}\sqrt{2\pi\sigma^2}}d\alpha \big)^{k}}{\big(\int_{-\frac{\bar{\tau}}{2}}^{\frac{\bar{\tau}}{2}}\frac{e^{-(o-a)^2/(2\sigma^2)}}{ \bar{\tau}\sqrt{2\pi\sigma^2}} d\alpha\big)^{k-1}} do.
\]

Then, to characterize the worst-case divergence, it is equivalent to considering the following

\begin{equation}\label{opt_h}
\max_{v} \sum_{j=1}^d h(v_j)
\text{~~ where ~~}
\sum_{j=1}^d v^2_j = \|v\|^2 \leq c^2.
\end{equation}

Due to the symmetric property of $h(\cdot)$, i.e., $h(v_j) = h(-v_j)$, we only need to consider the case where $v_j \geq 0$.
Define $g(x) = h(\sqrt{x})$, then solving Equation (<ref>) is equivalent to solving

\begin{equation}\label{opt_h2}
\max \sum_{j=1}^d g(\tilde{v}_j)
\text{~~ where ~~}
\forall j, \tilde{v}_j \in [0, c]
\text{~~ and ~~}
\sum_{j=1}^d \tilde{v}_j \leq c^2.
\end{equation}

With some simple calculation on the second-order derivative of $g(\tilde{v})$, $g(\tilde{v})$ is a convex function and therefore the maximal of $\sum_{j=1}^d g(\tilde{v}_j)$ over a simplex constraint must be achieved at the vertices.
Therefore, the maximum is achieved when $\tilde{v}$ is some one-hot vector, and we only need to consider the one-dimensional case by selecting $v = (c, 0, 0,...,0)$. As a straightforward proof of Corollary <ref>, when we further ensure the $l_{\infty}$ norm sensitivity, then $\tilde{v}$ is within the intersection of an $l_1$ ball of radius $c^2$ and a hyper cube, where the length of each side is $(c^2/p)$. Then, $\max_{v} \sum_{j=1}^d h(v_j)$ is achieved when $v  = (c/\sqrt{p}, ... ,c/\sqrt{p}, 0,...,0)$ whose Hamming weight is $p$. The rest proof is a straightforward application of Theorem <ref> to concert RDP results to $(\epsilon, \delta)$ DP to get the composition bound claimed.

When we apply Laplace mechanism to handle the case of $l_1$-norm sensitivity, i.e., $\sum_{j=1}^d |v_j| \leq c$, the corresponding $h(v_i)$ is still convex with respect to $v_i$. Therefore, the maximum of $\sum_{j=1}^d h(v_i)$ is still achieved when $v$ is a one-hot vector and we can reduce the multi-dimensional problem to the single dimension scenario.

§ PROOF OF THEOREM 3.3

Proof of Sketch:  Following the arguments of Theorem <ref>, still let $\mathcal{M}$ be the mechanism of one iteration of DP-SGD with ModelMix, and $w_1, w_2, ... , w_T$ be the outputs returned across $T$ iterations. To derive an $(\epsilon, \delta)$-DP bound of the cumulative privacy loss, it suffices to consider the high probability bound of the sum of point-wise privacy loss

$$\sup_{\mathcal{D},\mathcal{D}',\text{Aux}} \sum_{k=1}^T \epsilon_{\mathcal{D},\mathcal{D}',\mathsf{Aux}}(w_k|w_{[1:k-1]}),$$

for $\epsilon_{\mathcal{D},\mathcal{D}',\mathsf{Aux}}(w_k|w_{[1:k-1]}) = \log \frac{\mathbb{P}(\mathcal{M}(\mathcal{D}, \mathsf{Aux})= w_k|w_{[1:k-1]}) }{\mathbb{P}(\mathcal{M}(\mathcal{D'},\mathsf{Aux})= w_k|w_{[1:k-1]})}.$ In other words, to ensure $(\epsilon, \delta)$ DP of the cumulative privacy across $T$ iterations, it suffices to ensure

$$\Pr_{w_{[1:T]}}(\sup_{\mathcal{D},\mathcal{D}',\text{Aux}} \sum_{j=1}^T \epsilon_{\mathcal{D},\mathcal{D}',\mathsf{Aux}}(w_k|w_{[1:k-1]}) \leq \epsilon) \geq 1-\delta. $$

In our application, we may apply Bernstein-Azuma inequality <cit.> to derive a high probability bound of the sum of $\epsilon_{\mathcal{D},\mathcal{D}',\mathsf{Aux}}(w_k|w_{[1:k-1]})$ for $k=1,2,...,T$. Thus, we calculate the worst case expectation and the variance of $\epsilon_{\mathcal{D},\mathcal{D}',\mathsf{Aux}}(w_k|w_{[1:k-1]})$, respectively, and then apply Bernstein-Azuma inequality to  obtain the bound claimed. The details of calculation can be found below.

§.§ Gaussian mechanism case

Since we are interested in the asymptotic behavior of $\bar{\tau}$ and $q$ is a constant, it suffices to consider the case when $q=1$.
When $q$ = 1, we are essentially considering the full-batch GD and one can generalize the following analysis via the privacy amplification theorem by sampling <cit.>. Then, the two output distributions of (<ref>) from $\mathcal{D}$ and $\mathcal{D}'$ we aim to compare are equivalent to

$$ \mathbb{P}(o) = \int_{-\frac{\bar{\tau}}{2}}^{\frac{\bar{\tau}}{2}}\frac{1}{\bar{\tau}\sqrt{2\pi\sigma^2}} e^{-\frac{(o-\alpha-c)^2}{2\sigma^2}} d \alpha,$$

$$ \mathbb{P}'(o) = \int_{-\frac{\bar{\tau}}{2}}^{\frac{\bar{\tau}}{2}}\frac{1}{\bar{\tau}\sqrt{2\pi\sigma^2}} e^{-\frac{(o-\alpha)^2}{2\sigma^2}} d \alpha.$$

We define the pointwise $\epsilon(o)$ loss at $o$ as
$ \epsilon(o) = \log \frac{\mathbb{P}(o)}{\mathbb{P}'(o)}.$[Due to the symmetry of both uniform and Gaussian distributions, it is sufficient to only consider the case by defining $\epsilon(o) = \log \frac{\mathbb{P}(o)}{\mathbb{P}'(o)}$.] Regarding $\epsilon(o)$, we have the following observation. First,

\[
\mathbb{P}(o) = \frac{\Phi( \frac{o-c-\bar{\tau}/2}{\sigma})-\Phi(\frac{o-c+\bar{\tau}/2}{\sigma})}{\bar{\tau}},
\]

\[
\mathbb{P}'(o) = P(o + c) = \frac{\Phi(\frac{o-\bar{\tau}/2}{\sigma}) - \Phi( \frac{o+\bar{\tau}/2}{\sigma})}{\bar{\tau}},
\]

where $\Phi(t) = \int_{t}^{\infty} e^{-t^2/2}/\sqrt{2\pi}.$ Regarding $\Phi(t)$, we have the following folk lemma that for any $t >0$,

\begin{equation}\label{equ:errorF}
\frac{\sqrt{2}e^{-t^2/2}}{\sqrt{\pi}(t + \sqrt{t^2+4})} \leq  \Phi(t) \leq \frac{\sqrt{2}e^{-t^2/2}}{\sqrt{\pi}(t + \sqrt{t^2+8/\pi})}.
\end{equation}

For any $t < 0$, we can also bound $\Phi(t)$ since $\Phi(t) = 1 - \Phi(-t)$.

following observation. First, we note that the ratio

$$\gamma(\alpha, o,c)  = \frac{e^{-\frac{(o-\alpha)^2}{2\sigma^2}}}{e^{-\frac{(o-\alpha-c)^2}{2\sigma^2}}} = e^{\frac{-2(o-\alpha)c+c^2}{2\sigma^2}},$$

which increases as $\alpha$ increases, for any fixed $o$ and $c$. Therefore,

\begin{equation}
\begin{aligned}
\frac{-2(o+\bar{\tau}/{2})}{2\sigma^2}  = \log \gamma(-\bar{\tau}/2,&  o,c) \leq  \epsilon(o) \\
& \leq \log \gamma(\bar{\tau}/2, o,c)  = \frac{-2(o-\bar{\tau}/{2})}{2\sigma^2}
\end{aligned}
\end{equation}

that $\epsilon(o)$ increases when $\alpha \in (-\infty, c/2)$ and then decreases

\[
o \in \Big(s+ \frac{\log \tau}{\lambda}, \tau-\frac{\log \tau}{\lambda}\Big),
\]

To derive the composition bound of $(\epsilon ,\delta)$, we provide the following lemma to capture the mean and variance of $\epsilon(o)$, respectively.

When $\bar{\tau} > \max\{ 2\sigma\log \bar{\tau}+ 2c, e^{3c / \sigma}, 2\}$ and $c / \sigma = \Theta(1)$, for any $\rho$ such that

\[
\log\rho > {c\over \sigma} + \sqrt{2\log(\bar{\tau}\cdot \max\{c, 1\})},
\]

we have

\[
\begin{aligned}
\mathbb{E}_{o \sim P}~ \epsilon(o)
O({c\log^2\rho \over \bar{\tau}})
+ O({\sigma\log\rho\over \bar{\tau}}),
\end{aligned}
\]

and the variance $\mathcal{V}_0 = \text{Var}(\epsilon(o))$ satisfies

\[
\mathcal{V}_0
{1\over \bar{\tau}}\cdot O({c^4\log\rho\over \sigma^3} + {c^2\log^3\rho\over \sigma} + \sigma\log\rho).
\]

Let us consider the mean first.
Note that

\[
\begin{aligned}
\mathbb{E}_{o \sim P}~ \epsilon(o)
& =
\int P(o) \log({P(o)\over P'(o)})\\
& \le
\int_{P(o) > P'(o)} P(o) \log({P(o)\over P'(o)})\\
& =
\int_{c/2}^{+\infty} P(o) \log({P(o)\over P'(o)}).
\end{aligned}
\]

Therefore, we only consider the case when $o > c/2$.
Let $\rho$ denote a failure probability parameter that will be specified later.
We will divide the $[c/2, +\infty)$ integral into four parts
(1) $[c/2, \bar{\tau}/2 - \sigma\log\rho)$,
(2) $[\bar{\tau}/2 - \sigma\log\rho, \bar{\tau}/2 + 2c)$,
(3) $[\bar{\tau}/2 + 2c, \bar{\tau}/2 + \sigma\log\rho)$ and
(4) $[\bar{\tau}/2 + \sigma\log\rho, +\infty)$
and compute the integral of these four parts correspondingly.
When $o\in [c/2, \bar{\tau}/2 - \sigma\log\rho]$,

\[
\log{P(o)\over P'(o)}
= \log (1 + {P(o) - P'(o)\over P'(o)})
< {P(o) - P(o + c)\over P(o + c)}.
\]

By definition, we have

\[
P(o + c)
= {1\over \bar{\tau}}\cdot (1 - \Phi({-o + \bar{\tau}/ 2 \over \sigma}) - \Phi({o + \bar{\tau}/ 2 \over \sigma}))
> {1 - 2\cdot\Phi(\log\rho) \over \bar{\tau}},
\]

\[
\begin{aligned}
P(o) - P(o + c)
=~ &~
{1\over \bar{\tau}}\cdot
(\Phi({-o - c + \bar{\tau}/ 2 \over \sigma}) - \Phi({-o + \bar{\tau}/ 2 \over \sigma})) \\
&  + {1\over \bar{\tau}}\cdot
(\Phi({o + \bar{\tau}/ 2 \over \sigma}) - \Phi({o - c + \bar{\tau}/ 2 \over \sigma}))
\\
<~ & {1\over \bar{\tau}}\cdot \Phi({-o - c + \bar{\tau}/ 2 \over \sigma})
< {1\over \bar{\tau}}\cdot \Phi(\log \rho - {c \over \sigma}).
\end{aligned}
\]

Combining these with Equation (<ref>), we have that for any $o\in [c/2, \bar{\tau}/2 - \sigma\log\rho)$,

\[
\log{P(o)\over P'(o)}
< {\Phi(\log \rho- { c \over \sigma})
\over
1 - 2\cdot\Phi(\log \rho)}
< {e^{-(\log \rho - c/\sigma)^2 / 2}\over \log \rho - c/\sigma}.
\]

If we choose $\rho$ such that

\begin{equation}
\label{equ:rhoReq0}
\log \rho > {c\over \sigma} + \sqrt{2\log \bar{\tau}},
\end{equation}

\[
e^{-(\log \rho - c/\sigma)^2 / 2} < {1\over \bar{\tau}}
\implies
\log{P(o)\over P'(o)} < {1\over \bar{\tau}\cdot \sqrt{2\log \bar{\tau}}}.
\]

Thus, the integral of the first part $[c/2, \bar{\tau}/2 - \sigma\log\rho]$ is,

\[
\int_{c/2}^{\bar{\tau}/2 - \sigma\log\rho} P(o)\log{P(o)\over P'(o)} < {1\over \bar{\tau}\cdot \sqrt{2\log \bar{\tau}}}.
\]

Now, we consider the second part where $o\in [\tau/2 - \sigma\log\rho, \bar{\tau}/2 + 2c)$.
In this case, utilizing Equation (<ref>), we can show that

\[
\begin{aligned}
\log{P(o)\over P'(o)}
& <
\log{1 / \bar{\tau}\over P'(o)}
\le
\log{1\over \Phi({2c\over \sigma}) - \Phi({2c + \bar{\tau}/2 \over \sigma})}
O({c^2\over \sigma^2}).
\end{aligned}
\]

\[
\begin{aligned}
\int_{\tau/2 - \sigma\log\rho}^{\bar{\tau}/2 + 2c} P(o)\log{P(o)\over P'(o)} d o
& =
O({c^2\over \sigma^2}) \cdot {2c + \sigma\log\rho\over \bar{\tau}} \\
& =
O({c^2 \log\rho\over \sigma\bar{\tau}}).
\end{aligned}
\]

When $o \ge {\tau/2} + 2c$, we have

\[
\begin{aligned}
{P(o)\over P'(o)}
& =
{\Phi({o - c - \bar{\tau}/2 \over \sigma}) - \Phi({o - c + \bar{\tau}/2 \over \sigma})
\over
\Phi({o - \bar{\tau}/2 \over \sigma}) - \Phi({o + \bar{\tau}/2 \over \sigma})}
2 \cdot {\Phi({o - c - \bar{\tau}/2 \over \sigma})
\over
\Phi({o - \bar{\tau}/2 \over \sigma})} \\
%& <
%4 \cdot {o - \bar{\tau}/2\over o - c - \bar{\tau}/2}
%    \cdot e^{-{(o - c - \bar{\tau}/2)^2 / (2\sigma^2)} + {(o - \bar{\tau}/2)^2 / (2\sigma^2)}} \\
& <
4 \cdot e^{(2(o - \bar{\tau}/2)c - c^2) / (2\sigma^2)}.
\end{aligned}
\]

Here, the constant factor two in the first inequality is derived based on Equation (<ref>) and $\bar{\tau} > 1$.
When $\bar{\tau}$ is large, this factor will approximate 1.
We use 2 here as a simple relaxation.

\[
\log {P(o)\over P'(o)}
{2(o - \bar{\tau}/2)c - c^2 \over 2\sigma^2} + \log(4).
\]

For convenience, we will denote $z = \log(4) - c^2 / 2\sigma^2$ such that

\[
\log {P(o)\over P'(o)}
{(o - \bar{\tau}/2)c \over \sigma^2} + z.
\]

This allows us to compute the third and the fourth part as follows.
When $o\in [\bar{\tau}/2 + 2c, \bar{\tau}/2 + \sigma\log \rho)$,

\[
\begin{aligned}
& \int_{\bar{\tau}/2 + 2c}^{\bar{\tau}/2 + \sigma\log \rho} P(o)\log{P(o)\over P'(o)} d o \\
&  <
\int_{\bar{\tau}/2 + 2c}^{\bar{\tau}/2 + \sigma\log \rho} P(o)\cdot ({(o - \bar{\tau}/2)c\over \sigma^2} + z)  d o \\
& <
({c\log\rho\over \sigma} + z ) \cdot \int_{\bar{\tau}/2 + 2c}^{\bar{\tau}/2 + \sigma\log\rho} P(o) d o  \\
& <
({c\log\rho\over \sigma} + z )\cdot (\sigma\log\rho - 2c)\cdot P(\bar{\tau}/2 + 2c) \\
& =
O({c\log^2\rho \over \bar{\tau}}) + O({\sigma\log\rho\over \bar{\tau}}).
\end{aligned}
\]

Finally, when $o\in [\bar{\tau}/2 + \sigma\log \rho, +\infty)$,

\[
\begin{aligned}
\int_{{\bar{\tau}\over 2} + \sigma\log \rho}^{+\infty} P(o)\log{P(o)\over P'(o)} d o \\
& <
\int_{{\bar{\tau}\over 2} + \sigma\log \rho}^{+\infty} P(o)\cdot ({(o - \bar{\tau}/2)c\over \sigma^2} + z) do \\
& = O(
{1\over \bar{\tau}}\int_{{\bar{\tau}\over 2} + \sigma\log \rho}^{+\infty} e^{-(o - \bar{\tau}/2 - c)^2/(2\sigma^2)}\cdot {(o - \bar{\tau}/2)c\over \sigma^2}do ) \\
& = O(c
e^{-(\sigma\log \rho - c)^2/(2\sigma^2)}).
\end{aligned}
\]

In order for the integral of this part to be negligible, $\rho$ must be chosen such that

\begin{equation}
\label{equ:rhoReq1}
ce^{-(\sigma\log \rho - c)^2/(2\sigma^2)} < {1\over \bar{\tau}}
\implies
\log\rho > {c\over \sigma} + \sqrt{2\log(\bar{\tau}c)}.
\end{equation}

Combine this with the requirement in Equation (<ref>), we have the final requirement for $\rho$ that

\begin{equation}
\label{equ:rhoReq}
\log\rho > {c\over \sigma} + \sqrt{2\log(\bar{\tau}\cdot \max\{c, 1\})}.
\end{equation}

In conclusion, we can sum up the integral of the four parts to show that for any $\rho$ satisfying Equation (<ref>),

\[
\mathbb{E}_{o \sim P}~ \epsilon(o)
O({c^2 \log\rho\over \sigma\bar{\tau}})
+ O({c\log^2\rho \over \bar{\tau}})
+ O({\sigma\log\rho\over \bar{\tau}})
+ O({1\over \bar{\tau}}).
\]

Since $\log\rho > c / \sigma$, we can simplify the above equation to

\[
\mathbb{E}_{o \sim P}~ \epsilon(o)
=  O({c\log^2\rho \over \bar{\tau}})
+ O({\sigma\log\rho\over \bar{\tau}}).
\]

If we consider $c$ or $\sigma$ to be constant, then we can simply set $\rho = \Theta(\bar{\tau})$ to satisfy Equation (<ref>) and show that

\[
\mathbb{E}_{o \sim P}~ \epsilon(o)
= O({c\log^2\bar{\tau} \over \bar{\tau}})
+ O({\sigma\log\bar{\tau}\over \bar{\tau}}).
\]

Let us now consider the variance.

\[
\begin{aligned}
\mathcal{V}_0 (\epsilon(o))
= \mathbb{E}_{o \sim P} \epsilon^2(o)
= \int_{-\infty}^{c/2} P(o)\epsilon^2(o) do + \int_{c/2}^{+\infty} P(o)\epsilon^2(o) do
\end{aligned}
\]

Note that for any $o > c/2$,

\[
\begin{aligned}
\epsilon^2(o)
= \log^2{P(o)\over P'(o)}
= \log^2{P'({o})\over P(o)} \\
= \log^2{P'({-o})\over P(2c - o)}
= \log^2{P({c - o})\over P'(c - o)} \\
= \epsilon^2(c - o).
\end{aligned}
\]

However, due to property of the Gaussian function, for any $o > c/2$,

\[
P(o) = P(2c - o) > P(c-o).
\]

\[
\begin{aligned}
\int_{-\infty}^{c/2} P(o)\epsilon^2(o) do
& =
\int_{c/2}^{+\infty} P(c-o)\epsilon^2(c-o) do \\
& =
\int_{c/2}^{+\infty} P(c-o)\epsilon^2(o) do \\
& <
\int_{c/2}^{+\infty} P(o)\epsilon^2(o) do.
\end{aligned}
\]

Since we have already bounded $\epsilon(o)$ when computing the mean for any $o > c/2$,
we have

\[
\begin{aligned}
\mathcal{V}_0 (\epsilon(o))
& < 2 \int_{c/2}^{+\infty} P(o)\epsilon^2(o) do \\
& = {1\over \bar{\tau}}\cdot O({c^4\log\rho\over \sigma^3} + {c^2\log^3\rho\over \sigma} + \sigma\log\rho).
\end{aligned}
\]

The proof in cemma <ref> also implies that $\epsilon(o) < 1 / (\bar{\tau}\cdot \sqrt{2\log \bar{\tau}})$ in the first part $o\in [c/2, \bar{\tau}/2 - \sigma\log \rho)$.
Furthermore, the probability of $o$ being in the rest three parts is bounded by

\[
\Pr[o \ge \bar{\tau}/2 - \sigma\log \rho]
{2\sigma\log\rho \over \bar{\tau}}.
\]

We can choose $\rho$ to be the minimum that satisfies Equation (<ref>),
which implies that

\[
\Pr[o \ge \bar{\tau}/2 - \sigma\log \rho]
{2c + 2\sigma\sqrt{2\log (\bar{\tau} \cdot \max\{c, 1\})} \over \bar{\tau}}.
\]

\[
\begin{aligned}
\Pr\Big[\epsilon(o) \ge {1 \over \bar{\tau}\cdot \sqrt{2\log \bar{\tau}}}\Big]
& <
{2c + 2\sigma\sqrt{2\log (\bar{\tau} \cdot \max\{c, 1\})} \over \bar{\tau}}.
\end{aligned}
\]

For simplicity, let us define

\[
\epsilon_0 = {1 \over \bar{\tau}\cdot \sqrt{2\log \bar{\tau}}}
\text{~~and~~}
\delta_0 = {2c + 2\sigma\sqrt{2\log (\bar{\tau} \cdot \max\{c, 1\})} \over \bar{\tau}}.
\]

The proof in cemma <ref> also implies that $\epsilon(o)$ is small for the first three parts.
In other word, the probability that $\epsilon(o)$ is “large" is bounded by the probability that $o$ is in the fourth part $[\bar{\tau}/2 + \sigma\log \rho, +\infty)$.
Therefore, for any $\rho$ satisfying Equation (<ref>),

\[
\begin{aligned}
\Pr\Big[\epsilon(o) > \max({4c^2 \over \sigma^2}, {2c\log \rho \over \sigma})\Big]
& <
\Pr\Big[o \ge \bar{\tau}/2 + \sigma\log \rho\Big] \\
& <
{1\over \log \rho - c/\sigma}\cdot e^{-(\log \rho - c/\sigma)^2 / 2}.
\end{aligned}
\]

Since $\log \rho$ is at least $3c / \sigma$, we can rewrite the above inequality as

\[
\Pr\Big[\epsilon(o) > {2c\log \rho \over \sigma}\Big]
{3e^{-(2\log^2 \rho) / 9}\over 2\log \rho}.
\]

For simplicity, let us define

\[
\epsilon_0 = {2c\log \rho \over \sigma},
\text{~and~}
\delta_0 = {3e^{-(2\log^2 \rho) / 9}\over 2\log \rho}.
\]

The conditional expectation and variance satisfy

\[
\mathbb{E}_{o \sim P}~[\epsilon(o)~|~\epsilon(o) < \epsilon_0] = {\mathbb{E}_{o \sim P} \epsilon(o)\over 1 - \delta_0},
\]

\[
\mathbb{E}_{o \sim P}~[\epsilon^2(o)~|~\epsilon(o) < \epsilon_0] = {\mathcal{V}_0\over 1 - \delta_0}.
\]

Now, we may apply the Bernstein inequality <cit.> to give a concentration bound.
Let $o_1, o_2, ... ,o_T$ be the intermediate updates from Algorithm <ref> across $T$ iterations. Provided fresh randomness in each iteration, $\epsilon(o_k)$ forms an martingale. Conditional on that $\epsilon(o_k)$ that $|\epsilon(o_k)|\leq \epsilon_0$ for $k=1,2,...,T$, where we can apply Bernstein inequality, we have that with probability at least $1-(T\delta_0+ \tilde{\delta})$,

\begin{equation}
\begin{aligned}
\label{berstein_gau}
& \sum_{k=1}^T \epsilon(o_k)  = O\big( \frac{T\mathbb{E}[\epsilon(o)] + \sqrt{\log(1/\tilde{\delta})(T\mathcal{V}_0+\epsilon_0)} }{1-\delta_0} \big)\\
& =  \tilde{O}\big(\frac{T(c+\sigma)}{\bar{\tau}} + \sqrt{\log(1/\tilde{\delta})}(\sqrt{\frac{c^4/\sigma^3 + c^2/\sigma + \sigma}{\bar{\tau}}} + \frac{1}{\sqrt{\tau}})\big),
\end{aligned}
\end{equation}

for any $\tilde{\delta} \in (0,1)$. Finally, we scale the parameters back to the case of Algorithm <ref>, where the sensitivity is $c/B$ and $\bar{\tau}=\tau/\eta$, we have the bound claimed.

§.§ Laplace mechanism case

In the rest of the proof, we turn to consider the Laplace scenario with bounded $l_1$-norm sensitivity $c$. Now we assume that the noise $\Delta$ is a Laplace distribution of parameter $\lambda$, i.e., $\mathbb{P}(\Delta = o) = \frac{\lambda}{2} e^{-\lambda|o|}$. With a similar reasoning and by Remark <ref>, it is equivalent to considering the following distribution,

$$ o \sim \mathcal{U}[0,\bar{\tau}] + \text{cap}(\lambda).$$

With some calculation, we know its probability density function (pdf) can be expressed as follows,

\begin{equation*}
\begin{aligned}
\frac{e^{-\lambda |o|}(1-e^{-\lambda \bar{\tau}})}{2\bar{\tau}} & , & o \leq 0; \\
\frac{2-e^{-\lambda o} - e^{-\lambda(\bar{\tau} - o)}}{2\bar{\tau}}& , & o \in (0,\bar{\tau}); \\
\frac{e^{-\lambda |o-\bar{\tau}|}(1-e^{-\lambda \bar{\tau}})}{2\tau} & , & o \geq \bar{\tau}.
\end{aligned}
\right.
\end{equation*}

Similarly, we use $P'(o)$ to denote the pdf of  $\mathcal{U}[c,c+\tau] + \text{cap}(\lambda)$. We will use $\epsilon_0 = {\lambda c}$ in the following.

When $\bar{\tau} > \max\{ 2\log \bar{\tau}/\lambda+c,1\}$,

\[
\begin{aligned}\mathbb{E}_{o \sim P}~ \epsilon(o)
\epsilon_0(e^{\epsilon_0}-1)(\frac{1-e^{-\lambda \bar{\tau}}}{\bar{\tau}\lambda} + \frac{2\log \bar{\bar{\tau}}/\lambda  +c}{\bar{\tau}}) \\
+ \frac{e^{\epsilon_0}-1}{2(\bar{\tau}-1)}(e^{\frac{\epsilon_0}{2(\bar{\tau}-1)}}-1)( 1-\frac{1-e^{-\lambda \bar{\tau}}}{\bar{\tau}\lambda}).
\end{aligned}
\]

and the variance $\mathcal{V}_0 = \text{Var}(\epsilon(o))$ satisfies

\mathcal{V}_0  \leq
\epsilon^2_0(\frac{1-e^{-\lambda \bar{\tau}}}{\bar{\tau}\lambda} + \frac{2\log \bar{\tau}/\lambda +s}{\bar{\tau}}) + (\frac{e^{\epsilon_0}-1}{2(\bar{\tau}-1)})^2( 1-\frac{1-e^{-\lambda \bar{\tau}}}{\bar{\tau}\lambda}).

We split the calculation of $\mathbb{E}_{o \sim P} \epsilon(o)$ into two parts, where

\begin{equation*}
\begin{aligned}
\mathbb{E}_{o \sim P}~ \epsilon(o)
& = \int_{\mathbb{R}} P(o)\log \frac{P(o)}{P'(o)} ~ do \\
& \leq \int_{\mathbb{R}} P(o)\log \frac{P(o)}{P'(o)} ~ do + \int_{\mathbb{R}} P'(o)\log \frac{P'(o)}{P(o)} ~ do  \\
& = \int_{\mathbb{R}} \big [ P(o)(\log \frac{P(o)}{P'(o)} \\
& + \log \frac{P'(o)}{P(o)}) + (P'(o)-P(o))\log\frac{P'(o)}{P(o)} \big] do\\
& = \int_{\mathbb{R}}  (P'(o)-P(o))\log\frac{P'(o)}{P(o)}  do\\
& \leq \underline{\int_{\mathbb{R}/\mathcal{I}} \sup_{o \in \mathbb{R}/\mathcal{I}}|\log\frac{P'(o)}{P(o)}|\cdot|P'(o)-P(o)| do}_{(A)}\\
& \quad
+ \underline{\int_{\mathcal{I}} \sup_{o \in \mathcal{I}}|\log\frac{P'(o)}{P(o)}|\cdot|P'(o)-P(o)| do}_{(B)}.
\end{aligned}
\end{equation*}

Here, $\mathcal{I} = [c+\log \bar{\tau}/\lambda, \bar{\tau}-\log\bar{\tau}/\lambda]$.
We first handle (A), where it is easy to verify that ModelMix does not change the worst case $\epsilon$ bound and $ \sup_{o \in \mathbb{R}/\mathcal{I}}|\log\frac{P'(o)}{P(o)}| \leq \epsilon_0 = \lambda c$.
Thus, $|P'(o)-P(o)| \leq (e^{\epsilon_0}-1)P(o)$.
On the other hand, we can upper bound the probability $\Pr(o \in \mathbb{R}/\mathcal{I})$ as follows.
It can be computed that $\Pr(o \leq 0) = \Pr(o \geq \bar{\tau} ) = \frac{1-e^{-\lambda \bar{\tau}}}{2\bar{\tau}\lambda}.$
As for $o \in (0,\bar{\tau})$, $P(o)$ is indeed concentrated at its mean $\bar{\tau}/2$, and thus

\begin{equation}
\begin{aligned}
\Pr(o \in \mathbb{R}/\mathcal{I}) & \leq 2\cdot \frac{1-e^{-\lambda\bar{\tau}}}{2\bar{\tau}\lambda} + (1-\frac{|\mathcal{I}|}{\bar{\tau}})\cdot (1-2\frac{1-e^{-\lambda \bar{\tau}}}{2\bar{\tau}\lambda}) \\
& \leq \frac{1-e^{-\lambda \bar{\tau}}}{\bar{\tau}\lambda} + \frac{2\log \bar{\tau}/\lambda +c}{\bar{\tau}}.
\end{aligned}
\end{equation}

Therefore, we have $(A) \leq \epsilon_0(e^{\epsilon_0}-1)(\frac{1-e^{-\lambda \bar{\tau}}}{\bar{\tau}\lambda_0} + \frac{2\log \bar{\tau}/\lambda +c}{\bar{\tau}})$.

Now, we turn to bound (B). for sufficiently large

\[
\bar{\tau} > \max\Big\{ c+\frac{2\log \bar{\tau}}{\lambda}, 1\Big\},
\]

we have for $o \in \mathcal{I}$,

\begin{equation}
\label{middle}
\small{
\begin{aligned}
& = \big|\log\big(1 + \frac{e^{-\lambda(o-c)}+e^{-\lambda(\bar{\tau}+c-o)} - e^{-\lambda(o)}-e^{-\lambda(\bar{\tau}-o)}}{2-e^{-\lambda(o-c)}-e^{-\lambda(\bar{\tau}+c-o)}}\big)\big| \\
& = \big|\log\big(1 + \frac{ e^{-\lambda o}(e^{\lambda c}-1) + e^{-\lambda(\bar{\tau}-o)}(e^{-\lambda c}-1)}{2-e^{-\lambda(o-c)}-e^{-\lambda(\bar{\tau}+c-o)}}\big) \big| \\
& \leq \frac{(e^{\lambda c}-1)/\bar{\tau}}{2-2/\bar{\tau}} = \frac{e^{\epsilon_0}-1}{2\bar{\tau}-2}.
%\leq {\lambda s\over \tau}.
\end{aligned}
\end{equation}

\[
\sup_{o \in \mathcal{I}} \epsilon(o) \leq \frac{e^{\epsilon_0}-1}{2(\bar{\tau}-1)}
\]

$$\Pr(o \in \mathcal{I}) \leq 1-\Pr(o \in (-\infty,0)\cup(\bar{\tau}, +\infty)) = 1-\frac{1-e^{-\lambda \bar{\tau}}}{\tau\lambda},$$

we have

\[
\begin{aligned}\mathbb{E}_{o \sim P}~ \epsilon(o)
\epsilon_0(e^{\epsilon_0}-1)(\frac{1-e^{-\lambda \bar{\tau}}}{\bar{\tau}\lambda} + \frac{2\log \bar{\tau}/\lambda  +c}{\bar{\tau}}) \\
+ \frac{e^{\epsilon_0}-1}{2(\bar{\tau}-1)}(e^{\frac{\epsilon_0}{2(\bar{\tau}-1)}}-1)( 1-\frac{1-e^{-\lambda \bar{\tau}}}{\bar{\tau}\lambda}).
\end{aligned}
\]

In the following, we set out to characterize the variance of $\epsilon(o)$. It is noted that

\[
\begin{aligned}
\text{Var}(\epsilon(o))
& = \mathbb{E}\big[(\epsilon(o)-\mathbb{E}[\epsilon(o)])^2\big]
\leq \mathbb{E}[(\epsilon(o)-0)^2] \\
& = \mathbb{E}[\epsilon(o)^2\cdot (\bm{1}_{o \in \mathcal{I}} +\bm{1}_{o \in \mathbb{R}/\mathcal{I}}) ].
\end{aligned}
\]

We use the upper bound $|\epsilon(o)| \leq \frac{e^{\epsilon_0}-1}{2(\bar{\tau}-1)}$ when $o \in \mathcal{I}$ and for the rest we simply bound $|\epsilon(o)|\leq \epsilon_0$, and we have

\[
\begin{aligned}
\mathcal{V}_0 = \text{Var}(\epsilon(o))
\leq
\epsilon^2_0(\frac{1-e^{-\lambda \bar{\tau}}}{\bar{\tau}\lambda} + \frac{2\log \bar{\tau}/\lambda +s}{\bar{\tau}}) \\
& + (\frac{e^{\epsilon_0}-1}{2(\bar{\tau}-1)})^2( 1-\frac{1-e^{-\lambda \bar{\tau}}}{\bar{\tau}\lambda}).
\end{aligned}
\]

Still, applying Bernstein inequality, to ensure an $\delta$ failure probability, we may select

\[
\begin{aligned}
& = \sqrt{2} \cdot \sqrt{T\mathcal{V}_0\log(1/\delta)+\log(1/\delta)\epsilon^2_0/9} +\log(1/\delta)\epsilon^2_0/3 \\
& = \tilde{O}(\sqrt{T\log(1/\delta)} \cdot \big(\epsilon_0  \sqrt{\frac{1}{\bar{\tau}\lambda}+\frac{s}{\bar{\tau}}} + \frac{e^{\epsilon_0}-1}{\bar{\tau}} \big),
\end{aligned}
\]

and correspondingly,

\begin{equation}
\begin{aligned}
\epsilon & = \tilde{O}\Big( \epsilon_0\cdot {T\over \bar{\tau}}\cdot(e^{\epsilon_0}-1) + (\epsilon_0 + \frac{e^{\epsilon_0}-1}{\sqrt{\bar{\tau}}}) \cdot \sqrt{{T\over \bar{\tau}}\cdot \log(1/\delta)} \Big)\\
& = \tilde{O}(\frac{\eta T\epsilon_0(e^{\epsilon_0}-1)}{\tau} + \epsilon_0 \sqrt{\frac{\eta T \log(1/\delta)}{\tau}}).
\end{aligned}
\end{equation}

§ EXPERIMENTS SETUPS

Experiments in Section <ref> (The comparisons to <cit.> where we use DP-SGD to privately train a CNN on ScatterNetwork features.)
For CIFAR-10, we reproduce the results of <cit.> by choosing sampling rate $q=8192/50000$, clipping threshold $c=0.1$ with stepsize $\eta=4$ for $700$ iterations.
For FMNIST, to reproduce the results of <cit.>, we choose sampling rate $q=8192/60000$, clipping threshold $c=0.1$ with stepsize $\eta=4$ for $800$ iterations.
When we apply ModelMix, for both CIFAR10 and FMNIST, we adopt the same sampling rate as they use above.
But we use a larger clipping threshold $c=1$ with $\eta=0.4$ and run for the same number of iterations.
As for the selection of $\tau_k$, we select $\tau_k=0.05\eta$ for the first half epochs and $\tau_k=0.025\eta$ for the rest.

The CNN model that we use is the same as that applied in <cit.>. The parameters can be found in Table 8 and Table 9 in <cit.>.

Experiments in Section <ref> (The experiments on training Resnet20 with DP-SGD.)
For CIFAR-10, in the application of regular DP-SGD with per-sample clipping, we select the sampling rate $q=1500/50000$, clipping threshold $c=20$, $\eta=0.1$, and run for $3500$ iterations. When we apply ModelMix, we select the same $q=1500/50000$, $c=20$, but a larger $\eta=0.15$ and run for $3500$ iterations. We set $\tau_k=0.15\eta$ uniformly.

For SVHN, when we test regular DP-SGD, we select the sampling rate $q=1830/73260$, clipping threshold $c=20$, $\eta=0.1$, and run for $4000$ iterations. When we apply ModelMix, we select the same $q=1830/73260$, $c=20$, but a larger $\eta=0.15$ and run for $4000$ iterations. We set $\tau_k=0.15\eta$ uniformly.

Experiments in Section <ref>
(The experiments on DP-SGD applications with further access to public data.)
In the private transfer learning application, we reproduce the results of <cit.> on CIFAR-10 by applying DP-SGD to learn a linear model on the features transformed by pretrained SimCLRv2 network.
The parameters are selected as they suggest, where $q=1024/50000$, $c=0.1$ and $\eta=4$, and run for 2700 iterations.
When we incorporate ModelMix, we adopt the same $q=1024/50000$, $c=0.1$, but a larger stepsize $\eta=5$. $\tau$ is set as $0.0075\eta, 0.005\eta, 0.0025\eta$ for the first, second, and third 1/3 of the total 1500 iterations, respectively. Due to the privacy amplification of ModelMix, we use a smaller noise scaled by 0.85 compared to their selection. Consequently, we achieve $(\epsilon=0.64,\delta=10^{-5})$ with the same $92.7\%$ accuracy.

When comparing to the low-dimensional embedding method presented in <cit.>, we adopt $q=1000/50000$, $c_0=5$ for the embedded gradient and $c_1=2$ for the residual gradient, and $\eta=0.1$ for 10000 iterations, as suggested, to reproduce $73.2\%$ accuracy on CIFAR-10 with Resnet20 with a budget $(\epsilon=8,\delta=10^{-5})$ under Gaussian Mechanism. When we apply ModelMix, we select $q=2000/50000$, $\eta=0.2$, $p=25$, and with the same clipping thresholds as above and run for 5000 iterations. With the privacy amplification of ModelMix, we uniformly scale the noise added by $0.75$ and select $\tau=0.05$ uniformly for all iterations. We achieve an $74.2\%$ accuracy at a budget $(\epsilon =2.9, \delta=10^{-5})$. With the same setup,  we can also achieve an accuracy of $79.1\%$ at a budget $(\epsilon=6.1, \delta=10^{-5})$.