
Invertibility: \( Q + (\mu _P-\)

\(\phi _Q( = \phi _P(\)  and \(\phi _P( = \phi _Q(\)

\({\text{det}{\frac{d \ \text{reflect}(\mu _Q,\mu _P) }{d}}=1}\) 
[Proof of Lemma REF ]
The proofs follows immediately using the given assumptions.
Coupling Markov chains
We setup some notation to describe the coupling that alg:unlearn-noisy-m-a-sgd constructs.
The following discussion is for deletion of index \(n\) , but it can be verified that the arguments naturally extend to the insertion case.
We remind that \(\mu _{n,m}\)  denotes the distribution of sampling \(m\)  elements uniformly randomly from \([n]\) , and mini-batches \(b_j \sim \mu _{n,m}\) .
Furthermore, we will use \(\mathbf {b}_j = [b_1,b_2,\ldots ,b_j]\)  denote the set of indices upto \(j\) .
For dataset \(S\)  and mini-batch indices \(b\) , let the gradient \(\nabla \hat{F}_S(,_{b}) := \frac{1}{|b|}\sum _{j\in b} \nabla f(,_j)\) .
Define \(\tilde{}_{j+1} :=  \mathring{}_j - \eta \nabla F_S( \mathring{}_j,_{b_j})\) ,
\(\bar{}_{j+1} := \tilde{}_{j+1} -\eta \theta _j\)  and \(_{j+1} := (\bar{}_{j+1})\) . Note that \(\tilde{}_j\)  is also function of \(b_j\)  but this dependency is not highlighted for notational simplicity.
The iterates and the mini-batches \([(\bar{}_2,b_1),(\bar{}_3,b_2),\ldots , (\bar{}_{T+1},b_T)]\)  produced by Algorithm REF  is a sample from a \(T\) -step first order Markov Chain over an uncountable state space \(^d \times [n]^*\) .
We remark that \(\bar{}_1\)  is a constant initialization, and so isn't considered.
Let \(P\)  be the joint distribution over the \(T\)  iterates \(\times \)  mini-batches. The joint density of \(P\)  can be factored as,
\(\phi _P((\bar{}_2,b_1),(\bar{}_3,b_2),\ldots , (\bar{}_{T+1},b_T)) = \phi _{P}(\bar{}_2,b_1)\phi _{P}(\bar{}_3,b_2 | _2) \ldots \phi _{P}(\bar{}_{T+1},b_T|_{T},_{T-1})\) 

where \(\phi _P(\bar{}_2,b_1) = \bar{\phi }_P(\bar{}_2 | b_1)\mu _{n,m}(b_1)\)  and \(\bar{\phi }_P(\bar{}_2|b_1)\)  is the density of \((\tilde{}_2, \eta ^2\sigma ^2)\) .
Similarly, the conditionals \(\phi _P(\bar{}_j,b_{j-1}| _{j-1},_{j-2}) = \bar{\phi }_P(\bar{}_j | b_{j-1},_{j-1},_{j-2}) \mu _{n,m}(b_{j-1})\) .
Furthermore, let \(\tilde{P}\)  denote the marginal of \([\bar{}_2,\bar{}_2, \ldots , \bar{}_{T+1}]\) , the joint density of which can be factored as,
\(\tilde{\phi }_{\tilde{P}}(\bar{}_2, \bar{}_3, \ldots , \bar{}_{T+1}) = \phi _{\tilde{P}}(\bar{}_2)\phi _{\tilde{P}}(\bar{}_3 | _2) \ldots \phi _{\tilde{P}}(\bar{}_{T+1}| _{T},_{T-1})\) 

where \(\phi _{\tilde{P}}(\bar{}_2) = \mathbb {E}_{b_1}\phi _P(\bar{}_2,b_1)\) , and the conditional \(\phi _{\tilde{P}}(\bar{}_j| _{j-1}, _{j-2}) = \mathbb {E}_{b_{j-1}} \bar{\phi }_P(\bar{}_j, b_{j-1} | _{j-1},_{j-2})\) .
Finally, given a fixed mini-batch sequence \(\mathbf {b}={b_1,b_2,\ldots ,b_T}\) , let \(P_{\mathbf {b}}\)  denote the joint conditional distribution of \({\bar{}_2,\bar{}_3, \ldots , \bar{}_{T+1}}\)  given \(\mathbf {b}\) . In this case, \(P_{\mathbf {b}}\)  factorizes as:
\(\phi _{P_{\mathbf {b}}}(\bar{}_2, \bar{}_3, \ldots , \bar{}_{T+1}) = \phi _{P_{b_1}}(\bar{}_2)\phi _{ P_{b_2}}(\bar{}_3 | _2) \ldots \phi _{ P_{b_T}}(\bar{}_{T+1}| _{T},_{T-1})\) 

where \(\phi _{P_{b_1}}(\bar{}_2) = \bar{\phi }_P(_2 | b_1)\)  and \(\phi _{ P_{b_{j-1}}}(\bar{}_j | _{j-1},_{j-2}) = \bar{\phi }_P(\bar{}_j | b_{j-1}, _{j-1},_{j-2})\) .
We similarly have a Markov Chain to generate the iterates for dataset \(S^{\prime }\)  - call this joint distribution over iterates and mini-batches as \(Q\) , the marginals over iterates as \(\tilde{Q}\)  and for a given \(\mathbf {b} \sim {\mu _{n-1,m}}^{\otimes T}\) , the conditionals over the iterates as \(Q_{\mathbf {b}}\) .
<FIGURE>We now describe how the unlearning Algorithm REF  constructs a coupling between \(P\)  and \(Q\)  to generate \((\mathbf {\bar{w}}^{(1)},\mathbf {\bar{w}}^{(2)})\) . We first describe the coupling of mini-batch indices.
Sample \(\mathbf {b} \sim (\mu ^{m}_{n})^{\otimes T}\) , let \(\mathbf {b}^{(1)} = \mathbf {b}\) .
We now look at all \(b_j^{(1)} \in \mathbf {b}^{(1)}\) : if \(n \notin b_j^{(1)}\) , then let \(b_j^{(2)} = b_j^{(1)}\) , otherwise for each such \(b_j^{(1)}\) , we replace \(n\)  by randomly sampling an index from \([n] \backslash b_j^{(1)}\) , and call this \(b_j^{(2)}\) . We then define the ordered set \(\mathbf {b}^{(2)} = {b_j^{(2)}}_{j=1}^T\) .
From claim:sgddeletion, this is a valid coupling of mini-batch indices.
Sample \(\mathbf {\bar{w}} = [\bar{}_2,\bar{}_3,\ldots , \bar{}_{T+1}] \sim P_{\mathbf {b}^{(1)}}\) , which corresponds to training with Algorithm REF  on dataset \(S\) .
Set \(\mathbf {\bar{w}}^{(1)}:=\mathbf {\bar{w}}\) .
To generate \(\mathbf {\bar{w}}^{(2)}\) , we do rejection sampling steps at each iteration. At the first step, we sample \(u_1\sim \text{Unif}(0,1)\) , and check if \(u_1 \le \frac{\phi _{Q_{\mathbf {b}^{(2)}}}(\bar{}_2)}{\phi _{P_{\mathbf {b}^{(1)}}}(\bar{}_2)}\) . If the step succeeds, then we proceed to the second iteration, wherein we again do a step of rejection sampling with ratio of conditional densities and so on. However, if anyone of the rejection sampling step fails, lets say the \(t^{\text{th}}\)  step, then we do a reflection of iterate \(\bar{}_{t+1}\)  about the mid-point of the means of \(P_{\mathbf {b}^{(1)}}(\cdot |_{t},_{t-1})\)  and \(Q_{\mathbf {b}^{(2)}}(\cdot |_{t},_{t-1})\)  , which are \(\tilde{}^{(1)}_{t+1} = _{t} - \eta \nabla F_{S}(_{t},_{b_t})\)  and \(\tilde{}^{(2)}_{t+1} = _{t} - \eta \nabla F_{S^{\prime }}(_{t},_{b_t^{\prime }})\)  respectively.
Set \(\bar{}^{(2)}_{t+1} = \text{reflect}(\bar{}_{t+1},\tilde{}^{(2)}_{t+1},\tilde{}^{(1)}_{t+1})\) .
After the reflection, we continue training on dataset \(S^{\prime }\)  which corresponds to continue sampling from the \((t+1)^{\text{th}}\)  step of the Markov chain for \(Q_{\mathbf {b}^{(2)}}\)  conditioned on the \(t^{\text{th}}\)  sample being \(\bar{}^{(2)}_{t}\) .
This generates the random variables \(\mathbf {\bar{w}}^{(1)}\)  and \(\mathbf {\bar{w}}^{(2)}\) .
We now show that this is indeed a coupling.
For any measurable set \(E \subseteq ^{dT}\) , \(¶{\mathbf {\bar{w}}^{(2)} \in E} = \tilde{Q}(E)\) 

[Proof of lem:m-sgdcoupling]
We will first show that \(¶{\mathbf {\bar{w}}^{(2)}\in E \vert (\mathbf {b}^{(1)}, \mathbf {b}^{(2)})} = Q_{\mathbf {b}^{(2)}}(E)\) . The proof is based on induction on the length of the Markov chain \(T\) .
Define \(\mathbf {\bar{w}}_T^{(2)} := {\bar{}_2^{(2)},\cdots , \bar{}_{T-1}^{(2)}}\) .
The key to the proof is the observation that the marginals \({P_{b_1^{(1)}}}(\cdot )\)  and \({Q_{b_1^{(2)}}}(\cdot )\)  are Gaussian \((\tilde{}_2^{(1)},\eta ^2\sigma ^2)\)  and \((\tilde{}_2^{(2)},\eta ^2\sigma ^2)\)  respectively, and the conditionals \({P_{b_j^{(1)}}}(\cdot | _{j})\)  and \({Q_{b_j^{(2)}}(\cdot | _{j})}\)  are also Gaussian \((\tilde{}_{j+1}^{(1)},\eta ^2\sigma ^2)\)  and \((\tilde{}_{j+1}^{(2)},\eta ^2\sigma ^2)\) .
For \(T=1\) , we only care about the marginals \({P_{b_1^{(1)}}}(\cdot )\)  and \({Q_{b_1^{(2)}}}(\cdot )\) , which as argued before, are normally distributed.
From Lemma REF , we have established that the reflection map satisfies the conditions in Lemma REF .
Combining these,
we have that the base case \(T=1\)  follows from the reflection coupling result stated as Lemma REF .
We proceed to the induction step.
There are two cases, depending on whether we do a rejection sampling in the \(T^{th}\)  step or not: we call these "rej-sample" and "no-rej-sample" respectively. If we do a rejection sampling, we further have two cases (1a). accept: either all rejection samplings, including the one in the \(T^{\text{th}}\)  step are accepts, (1b). reflect: all rejection samplings, except the one in the \(T^{\text{th}}\)  step are accepts, and in the \(T^{\text{th}}\)  step, we reflect. Finally, if we don't do a rejection sampling step, we have the third case (2). reject: some rejection sampling prior to \(T\)  results in reject; in this case, the \(T^{\text{th}}\)  sample \(\bar{}_{T+1}^{(2)} \sim Q_{b_T^{(2)}}(\cdot | _{T}^{(2)}, _{T-1}^{(2)})\) .
Cases (1) and (2) partition the whole event space for \(T\)  draws, whereas cases (1a) and cases (1b) partitions the space of the \(T^{\text{th}}\)  draw, conditioned on the first event. Also note that case (1) vs (2) distinction is measurable w.r.t. the natural filtration generated by the Markov chain upto \(T-1\)  draws.
Note that conditioned on the events "rej sample" as well as \(\mathbf {w}^{(2)}_{T-1}\) , the last step is just a one-step reflection coupling method. To elaborate, the conditionals \(Q_{b_T^{(2)}} (\cdot | _{T}^{(2)},_{T-1}^{(2)})\)  and \(P_{b_T^{(2)}} (\cdot | _{T}^{(2)},_{T-1}^{(2)})\)  used in the \(T^{\text{th}}\)  rejection sampling are Gaussians, which along with the reflection map satisfies properties of Lemma REF , as in the base case.
Let \(E_{last} = {y | \exists \mathbf {x}: (\mathbf {x},y) \in E}\)  be the projection of \(E\)  on the last co-ordinate and \(E_{y} = { \mathbf {x}| (\mathbf {x},y) \in E}\) .
According to lem:reflection, the conditional distribution of \(_{T+1}^{(2)}\)  is \(Q_{b_T^{(2)}} (\cdot | _{T}^{(2)},_{T-1}^{(2)})\) :
\(&¶{\bar{}_{T+1}^{(2)}\in E_{last} \vert \text{rej-sample}, { \mathbf {\bar{w}}_{T-1}^{(2)}}, (\mathbf {b}^{(1)},\mathbf {b}^{(2)})} \\&= ¶{\bar{}_{T+1}^{(2)}\in E_{last}, \text{accept} \vert \text{rej-sample}, { \mathbf {\bar{w}}_{T-1}^{(2)}}, (\mathbf {b}^{(1)},\mathbf {b}^{(2)})} \\&\qquad +¶{\bar{}_{T+1}^{(2)}\in E_{last}, \text{reflect} \vert \text{rej-sample}, { \mathbf {\bar{w}}_{T-1}^{(2)}}, (\mathbf {b}^{(1)},\mathbf {b}^{(2)})} \\&= \Large \int _{^{d}} \mathbb {1}{\bar{}_{T+1}^{(2)}\in E_{last}}\phi _{Q_{b_T^{(2)}}}(\bar{}_{T+1}^{(2)}| _{T}^{(2)}, _{T-1}^{(2)})d\bar{}_{T+1}^{(2)}\) 
For the "no-rej-sample" case, we have:
\(&¶{\bar{}_{T+1}^{(2)}\in E_{last}\vert \text{no-rej-sample}, { \mathbf {\bar{w}}_{T-1}^{(2)}}, (\mathbf {b}^{(1)},\mathbf {b}^{(2)})} =\underset{\bar{}_{T+1}^{(2)}}{\mathbb {E}} \left[ \mathbb {1}{\bar{}_{T+1}^{(2)}\in E_{last}}\right]\\& = \Large \int _{^{d}} \mathbb {1}{\bar{}_{T+1}^{(2)}\in E_{last}} \phi _{Q_{b_T^{(2)}}}(\bar{}_{T+1}^{(2)}|_{T}^{(2)},_{T-1}^{(2)}) d\bar{}_{T+1}^{(2)}\\\) 
We will now combine the two cases. Let \(\phi ^{(r)}(\cdot )\)  and \(\phi ^{(nr)}(\cdot )\)  denote the densities of \(\mathbf {\bar{w}}_{T-1}^{(2)}\)  under the "rej-sample" and "no-rej-sample" events respectively.
\(&¶{\mathbf {\bar{w}}^{(2)} \in E | (\mathbf {b}^{(1)},\mathbf {b}^{(2)})}\\&= ¶{\bar{}_{T+1}^{(2)}\in E_{last} \, \Big \vert \, \text{rej-sample}, { \mathbf {\bar{w}}_{T-1}^{(2)}}, (\mathbf {b}^{(1)},\mathbf {b}^{(2)})} ¶{\mathbf {\bar{w}}_{T-1}^{(2)} \in E_{\bar{}_{T+1}^{(2)}},\text{rej-sample}}\\&+ ¶{\bar{}_{T+1}^{(2)}\in E_{last} \, \Big \vert \, \text{no-rej-sample}, { \mathbf {\bar{w}}_{T-1}^{(2)}}, (\mathbf {b}^{(1)},\mathbf {b}^{(2)})} ¶{\mathbf {\bar{w}}_{T-1}^{(2)} \in E_{\bar{}_{T+1}^{(2)}},\text{no-rej-sample}}\\& = \int _{^{dT}}\mathbb {1}{\bar{}_{T+1}^{(2)}\in E_{last}}\mathbb {1}{\mathbf {\bar{w}}_{T-1}^{(2)} \in E_{\bar{}_{T+1}^{(2)}}} \\&\cdot {\mathbb {1}_\text{rej-sample}{\mathbf {\bar{w}}_{T-1}^{(2)}} \phi ^{(r)}(\mathbf {\bar{w}}_{T-1}^{(2)}) + \mathbb {1}_{\text{no-rej-sample}}{\mathbf {\bar{w}}_{T-1}^{(2)}}\phi ^{(nr)}(\mathbf {\bar{w}}_{T-1}^{(2)})}\phi _{Q_{{b_T^{(2)}}}}(\bar{}_{T+1}^{(2)}|_{T}^{(2)},_{T-1}^{(2)}) d \mathbf {\bar{w}}_{T}^{(2)}\\& = \int _{^{dT}}\mathbb {1}{\mathbf {\bar{w}}^{(2)} \in E}\phi _{Q_{\mathbf {b}_{T-1}^{(2)}}}(\mathbf {\bar{w}}_{T-1}^{(2)})\phi _{Q_{{b^{(2)}_T}}}(\bar{}_{T+1}^{(2)}| _{T}^{(2)},_{T-1}^{(2)}) d \mathbf {\bar{w}}_{T}^{(2)}\\& = \int _{^{dT}}\mathbb {1}{\mathbf {\bar{w}}^{(2)} \in E}\phi _{Q_{\mathbf {b}^{(2)}}}(\mathbf {w}_{T}^{(2)})d \mathbf {\bar{w}}_{T}^{(2)} = Q_{\mathbf {b}^{(2)}}( E)\) 

where the third equality uses the induction hypothesis that \(\mathbf {\bar{w}}_{T-1}^{(2)}\) , conditioned on \(\mathbf {b}^{(1)}\)  and \(\mathbf {b}^{(2)}\) , is distributed as \(Q_{\mathbf {b}_{T-1}^{(2)}}\) .
Finally, we integrate with respect to the coupling generating \( (\mathbf {b}^{(1)},\mathbf {b}^{(2)})\) ; we get
\(¶{\mathbf {\bar{w}}^{(2)} \in E} &= \sum _{(\mathbf {b}^{(1)},\mathbf {b}^{(2)})}¶{\mathbf {\bar{w}}^{(2)} \in E | (\mathbf {b}^{(1)},\mathbf {b}^{(2)})} ¶{\mathbf {b}^{(1)},\mathbf {b}^{(2)}}\\& =\sum _{(\mathbf {b}^{(1)},\mathbf {b}^{(2)})}Q_{\mathbf {b}^{(2)}}( E) ¶{\mathbf {b}^{(1)},\mathbf {b}^{(2)}} = \sum _{\mathbf {b}^{(2)}}Q_{\mathbf {b}^{(2)}}( E) ¶{\mathbf {b}^{(2)}} = \tilde{Q}(E)\) 
This completes the proof.

We now show that not only the marginals over the iterates, but the entire state maintained by the algorithm, which includes the mini-batching indices is transported.


For any measurable event in \({^d \times [n]^m}^{\otimes T}\) , we have
\(¶{(\mathbf {\bar{w}}^{(2)}, \mathbf {b}^{(2)}) \in E} = Q(E)\) 


[Proof of lem:noisy-m-sgd-transport]
We first decompose the event \(E \subseteq {^d \times [n]^m}^{\otimes T}\)  as two events, \(E=E_1 \times E_2\)  where \(E_1 \subseteq ^{dT}\)  and \(E_2 \subseteq {[n]^m}^T\) . We have
\(¶{(\mathbf {\mathbf {\bar{w}}}^{(2)}, \mathbf {b}^{(2)}) \in E} &=\mathbb {E}_{\mathbf {b}^{(1)}}¶{\mathbf {\mathbf {\bar{w}}}^{(2)} \in E_1 |\mathbf {b}^{(1)}, \mathbf {b}^{(2)} } ¶{\mathbf {b}^{(2)} \in E_2} \\& = \mathbb {E}_{\mathbf {b}^{(1)}} \tilde{Q}_{\mathbf {b}^{(2)}}(E_1)\mu _{n,m}^{\text{del}\otimes T}(E_2) \\& = \mathbb {E}_{\mathbf {b}^{(1)}} \tilde{Q}_{\mathbf {b}^{(2)}}(E_1)\mu _{n-1,m}^{\otimes T}(E_2) \\& = \tilde{Q}_{\mathbf {b}^{(2)}}(E_1)\mu _{n-1,m}^{\otimes T}(E_2) = Q(E)\) 
where the second and third equality follows from lem:m-sgdcoupling and claim:sgddeletion, and the final equality follows from the definition of event \(E\)  and probability distribution \(Q\) .
We now lower bound the probability of accepting at all rejection sampling steps.
Let “accept" be the event in which all rejection sampling result in accepts so there is no reflection or recompute. The probability of accept is lower bounded as,
\(¶{\text{accept}} \ge 1 - \frac{\sqrt{T}\rho }{8}\) 

[Proof of lem:noisy-m-sgd-accept]
We evaluate the probability that all rejection sampling steps result in accepts. We first do it conditioned on \(\mathbf {b}^{(1)},\mathbf {b}^{(2)}\)
\(¶{\text{accept}| \mathbf {b}^{(1)},\mathbf {b}^{(2)}} &= \underset{\bar{}_T^{(1)}, \bar{}_{T+1}^{(2)}, {u_j}_{j=1}^T}{\mathbb {E}} \left[ \prod _{j=1}^T\mathbb {1}_{\text{accept}}(u_j)\right]\\& = \underset{\mathbf {\bar{w}}_T^{(1)}, \mathbf {\bar{w}}_{T+1}^{(2)}}{\mathbb {E}} \left[ \prod _{j=1}^T ¶{u_j \le \frac{\phi _{Q_{b^{(2)}_j}} (\bar{}_{j+1}^{(1)} | _{j}^{(1)}, _{j-1}^{(1)}) }{\phi _{P_{b^{(1)}_j}} (\bar{}_{j+1}^{(1)} | _{j}^{(1)}, _{j-1}^{(1)}) } \Big \vert \mathbf {\bar{w}}_T^{(1)}, \mathbf {\bar{w}}_T^{(2)}}\right]\\& = \int _{^{dT}} \prod _{j=1}^T \min {\phi _{P_{b^{(1)}_j}} (\bar{}_{j+1}^{(1)} | _{j}^{(1)}) _{j-1}^{(1)}) , \phi _{Q_{b^{(2)}_j}} (\bar{}_{j+1}^{(1)} | _{j}^{(1)}), _{j-1}^{(1)})} d \mathbf {\bar{w}}^{(1)}_{T} \\& = \prod _{j=1}^T {\int _{^{d}} \min {\phi _{P_{b^{(1)}_j}} (\bar{}_{j+1}^{(1)} | _{j}^{(1)}), _{j-1}^{(1)}) , \phi _{Q_{b^{(2)}_j}} (\bar{}_{j+1}^{(1)} | _{j}^{(1)}), _{j-1}^{(1)})} d \bar{}^{(1)}_{j+1} }\\\) 
The term
\(&\int _{^{d}} \min {\phi _{P_{b^{(1)}_j}} (\bar{}_{j+1}^{(1)} | _{j}^{(1)}), _{j-1}^{(1)}) , \phi _{Q_{b^{(2)}_j}} (\bar{}_{j+1}^{(1)} | _{j}^{(1)}), _{j-1}^{(1)})} d \bar{}^{(1)}_{j+1} \\&= 1 -\text{TV}_{\mathbf {\bar{w}}_{j}^{(1)}, \mathbf {b}^{(1)},\mathbf {b}^{(2)} }(P_{b_j^{(1)}},Q_{b_j^{(2)}})\) 

where the notation \(\text{TV}_{x}(\cdot , \cdot )\)  denotes the conditional TV between the arguments, conditioned on the subscript.
Let \(\gamma _j= \frac{\Delta (b_j^{(1)},b_j^{(2)})}{2}\)  i.e. the number of elements differing in \(b_j^{(1)}\)  and \(b_j^{(2)}\)  . Note that if \(\gamma _j = 0\) , then \(b_j^{(1)}= b_j^{(2)}\)  and \(P_{b_j^{(1)}} = Q_{b_j^{(2)}}\) , and hence \(\text{TV}_{\bar{}_{j-1}^{(1)}, \mathbf {b}^{(1)},\mathbf {b}^{(2)} }(P_{b_j^{(1)}},Q_{b_j^{(2)}}) = 0\) . In the other case, \(\gamma _j = 1\) , which corresponds to the case when the deleted point was used in the \(j^{\text{th}}\)  mini-batch. In this case, the means of \(P_{b_j^{(1)}}\)  are \(Q_{b_j^{(2)}}\)  at separated by at most \(\frac{2G\eta }{m}\)  - this follows as in the proof of Proposition REF . In particular, fixing previous iterates and \(\bar{}_{j-1}^{(1)}\)  and mini-batch indices \(\mathbf {b}^{(1)},\mathbf {b}^{(2)}\) , using the fact that gradients are in norm bounded by \(G\) , \(P_{b_j^{(1)}}\)  and \( Q_{b_j^{(2)}}\)  are Gaussians with variance \(\eta ^2\sigma ^2\) 
and means separated by either \(\frac{2G\eta \gamma _j}{m}\)  or 0, depending on \(\gamma _j\) . Combining the two cases, and using TV between Gaussians formula [1]}, we have \(1- \text{TV}_{\bar{}_{j-1}^{(1)}, \mathbf {b}^{(1)},\mathbf {b}^{(2)} }(P_{b_j^{(1)}},Q_{b_j^{(2)}}) \ge { 1- \frac{G\eta }{\eta \sigma m}}^{\gamma _j} = { 1- \frac{G}{\sigma m}}^{\gamma _j}\) . We therefore get \(\int _{^{d}} \min {\phi _{P_{b^{(1)}_j}} (\bar{}_{j+1}^{(1)} | _{j}^{(1)}), _{j-1}^{(1)}) , \phi _{Q_{b^{(2)}_j}} (\bar{}_{j+1}^{(1)} | _{j}^{(1)}), _{j-1}^{(1)})} d \bar{}^{(1)}_{j+1} \ge { 1- \frac{G}{\sigma m}}^{\gamma _j}\) . Plugging this in the conditional probability of accept expression, we get
\(¶{\text{accept}| \mathbf {b}^{(1)},\mathbf {b}^{(2)}} &\ge \prod _{j=1}^T { 1- \frac{G}{\sigma m}}^{\gamma _j} = {1- \frac{G}{\sigma m}}^{\sum _{j=1}^T \gamma _j} \ge 1 - \frac{G\sum _{j=1}^T \gamma _j}{\sigma m}\) 
We now integrate with respect to \(\mathbf {b}^{(1)},\mathbf {b}^{(2)}\) . Note that \(\sum _{j=1}^T \gamma _j\)  is the number of mini-batches which contain the deleted point. Since in each mini-batch, \(m\)  points are selected uniformly randomly from \(n\) , \(\mathbb {E}_{\mathbf {b}^{(1)},\mathbf {b}^{(2)}}{\gamma _j} = \frac{m}{n}\) , which gives us \(\mathbb {E}_{\mathbf {b}^{(1)},\mathbf {b}^{(2)}}\sum _{j=1}^T \gamma _j = \frac{Tm}{n}\) . Hence,
\(¶{\text{accept}} \ge 1 - \frac{G \mathbb {E}_{\mathbf {b}^{(1)},\mathbf {b}^{(2)}}\sum _{j=1}^T \gamma _j}{\sigma m} = 1- \frac{GT}{\sigma n} = 1 -\frac{\sqrt{T}\rho }{8}\) 

where the last equality follows from plugging in \(\sigma = \frac{8\sqrt{T}G}{n \rho }\)  as in prop:upper-bound-noisy-m-a-sgd.
We are now ready to prove prop:unlearn-noisy-m-A-SGD.
[Proof of prop:unlearn-noisy-m-A-SGD]
We need to show that upon deletion and insertion, the probability distribution of the entire state maintained by the algorithm, which is all iterates as well as mini-batches indices is transported - this, for one deletion, follows from Lemma REF  (which is for unprojected iterates), together with the fact that projection is a deterministic operation.
Moreover, as before, the above argument also holds for insertion and generalizes arbitrary edit requests.
We now proceed to bound the probability to recompute.
This follows directly combining lem:noisy-m-sgd-accept, prop:upper-bound-noisy-m-a-sgd and Remark REF .
From Remark REF , upon \(k\)  edits, the total variation distance is at most \(k\)  times total variation distance upon 1 edit.
Since the algorithm is \(\rho \) -TV stable (Proposition REF ), and the assumption that the number of samples are between \(n/2\)  and \(2n\) , the total variation distance upon \(k\)  edits is at most \(2k\rho \) .
Hence, using lem:noisy-m-sgd-accept, we have that the probability to recompute is probability of “reject" is at most \(\frac{k\rho \sqrt{T}}{4}\) .

Runtime and space complexity
In this section, we discuss, in detail, the learning and unlearning runtime of the algorithms, as well as their space complexity.
Learning runtime
In this work, we did not aim to carefully optimize the runtime for training/learning algorithm, as long as the algorithm achieves the rate in thm:main-upper-bound.
However, we briefly discuss the runtime of each algorithm, and highlight easy improvements, where possible. Algorithm REF  requires \(mT\)  = \(\rho n\)  stochastic gradient computations.
On the other hand, for Algorithm REF , if \(m\le \sqrt{\frac{LD \rho n}{G \sqrt{d}}}\) , it requires \(mT = m{\frac{\rho n}{md}} = \frac{(\rho n)^2}{d}\)  stochastic gradient computations - setting larger \(m\)  only hurts the total runtime, without any advantage. Note that total stochastic gradient descent computations of noisy-m-SGD (i.e. without acceleration, see sec:noisy-m-sgd) is also \(\frac{(\rho n)^2}{d}\) ; however, the key advantage of acceleration is that it allows setting larger mini-batch sizes: \(T^3\)  as opposed to \(T\) , which leads to smaller number of iterations: \(\sqrt{\frac{\rho n}{\sqrt{d}}}\)  as opposed to \(\frac{\rho n}{\sqrt{d}}\)  and hence a smaller probability of recompute.
From [2]}, we know that mini-batch SGD (with or without acceleration) is optimal for smooth convex composite/ERM optimization in the low accuracy regime: when accuracy \(\alpha \gtrsim \frac{1}{\sqrt{n}}\) . In this regime, an algorithm makes at least \(\Omega {\frac{1}{\alpha ^2}}\)  calls to a stochastic gradient oracle.
It can be then verified that for our accuracy, our algorithms make the optimal number of oracle calls.
For Algorithm REF , as discussed,
faster algorithms lead to better unlearning times.
It is natural to ask what happens if we additionally introduce variance reduction techniques on top of acceleration to yield even faster runtimes.
In particular, what if we use Katyusha [3]}, which has optimal runtime in terms of stochastic gradient computations.
We argue that even though it improves the runtime of the learning algorithm, it does not yield improvement for unlearning beyond what we have with acceleration.
From Corollary 5.8 in [3]}, setting largest allowed \(m=\sqrt{n}\) , we get that \(Tm = O{\frac{\sqrt{n}}{\sqrt{\epsilon }}}\)  – in our case, \(\epsilon = \frac{\sqrt{d}}{\rho n}\) , which yields \(Tm = O{\sqrt{n}\sqrt{\frac{\rho n}{\sqrt{d}}}}\)  stochastic gradient computations. Note that this is smaller than that of noisy-m-A-SGD (unless \(\rho \)  is very small), however \(T = \sqrt{\frac{\rho n}{\sqrt{d}}}\)  – same as that of noisy-m-A-SGD, and hence yields no improvement in unlearning time.
However, note that using Katyusha would give us optimal oracle complexity even in the high accuracy regime.

Unlearning runtime
We now look at how much compute it takes for Algorithm REF  and REF  to handle the edit requests.
We first give a general result, which holds for any TV-stable algorithm with the unlearning algorithm being the one which constructs a coupling with acceptance probability at least \(1-\rho \) .
We give in-expectation bounds
on the number of times verification fails or a full or partial recompute is triggered.
For a coupling based unlearning algorithm with acceptance probability at least \(1-\rho \) , for \(k\)  edit requests, the expected number of times recompute is triggered is at most \(4k \rho \) .

[Proof of Proposition REF ]
We first setup some notation.
In the general setup, for \(k\)  edit requests, let \(s\)  be the number of times a recompute is triggered.
Let \({Z_1,Z_2,\ldots , Z_s}\)  be a set of random variables, where each \(Z_i\)  denotes how many edit requests the \(i^{\text{th}}\)  recompute can handle.
To elaborate, \(Z_i\)  takes value \(j\) , if upon \(j\)  edit requests, a recompute is triggered.
The \(Z_i\) 's comprises to the randomness used in the algorithm like mini-batching indices or Gaussian noise, as well as the randomness used for rejection sampling.
It is important to note that \(Z_i^{\prime }s\)  are not necessarily independent.
In particular, in Algorithm REF , we reuse the Gaussian noise upto the iteration in which rejection sampling fails, and only use fresh/independent Gaussian noise in the later steps.
However, note that we have exact unlearning, and the output at each step is \(\rho \) -TV stable (w.r.t. all the randomness used). Hence, since the above description of the distribution of \(Z_i\) 's depend only on the TV stability parameter,
it follows that \(Z_i^{\prime }s\)  are (marginally) identical.
We now use the fact that the unlearning algorithm constructs a coupling with acceptance probability at least \(1-\rho \)  to describe the probability distribution of \(Z_i\) .
We have that upon one edit request, the probability that a recompute is triggered is at most \(\rho \) .
This means that \(Z_i < 2\)  with probability \(\le \rho \) .
Using Remark REF , this generalizes as \(Z_i < j\)  with probability at most \((j-1)\rho \) .
Note that in our setup, we observe at most \(k\)  requests, so \(Z_i\)  taking values larger than \(k\)  is not meaningful. However if \(k\rho < 1\) , it means that probability that \(Z_i\)  takes values smaller than \(k\)  is less than 1, and therefore there is a positive probability of \(Z_i\)  being larger than \(k\) .
To remedy this, we define another random variable \(X_i\) 's which takes values in the set \({1,2,\ldots , k}\) .
Furthermore, for any \(i\) , \(¶{X_i = j} = ¶{Z_i = j}\)  for \(j\le k\) , but \(¶{Z_i = j} = \sum _{l=k}^\infty ¶{X_i=l}\) . By construction, this ensures that \(1\le Z_i\le k\) , when we observe at most \(k\)  requests.
We want upper bounds on \(s\)  conditioned on the fact that \(k\)  requests are addressed i.e. \(X_1+X_2 \ldots , X_s \ge k\) .
For this we write \(s\)  as \(s:=\min _q{2k \ge X_1+X_2 \ldots , X_q \ge k}\) .
The first inequality holds trivially since we ensured that \(X_i\le k\) .
It is easy to see that \(s\)  is a stopping time with respect to the filtered probability space of the stochastic process \({X_i}_{i\in }\) .
Furthermore, since \(X_i\) 's are identical, we can apply Wald's equation to get,
\(&2k \ge \mathbb {E}[X_1+X_2 + \ldots X_s] = \mathbb {E}[s] \mathbb {E}[X_1] = \mathbb {E}[s] \sum _{j=1}^{k} {X_i \ge j} =\mathbb {E}[s] \sum _{j=1}^{\infty } {Z_i \ge j}\\& \ge \mathbb {E}[s] \sum _{j=1}^{1+1/\rho } {Z_i \ge j} = \mathbb {E}[s] \sum _{j=1}^{1+1/\rho }{1- {Z_i < j}} \ge \mathbb {E}[s]{\frac{1}{\rho } - \sum _{j=1}^{1+1/\rho }(j-1)\rho } \\&\ge \mathbb {E}[s]{\frac{1}{\rho } - \int _{j=0}^{1/\rho }j\rho dj} \ge \mathbb {E}[s]{\frac{1}{\rho } - \frac{1}{2\rho ^2}\rho }\ge \frac{\mathbb {E}[s]}{2\rho }\) 

This gives us that \(\mathbb {E}[s] \le 4k\rho \) .
Next, we look at the runtimes for Algorithm REF  and Algorithm REF  to handle one deletion or insertion request. For this, we look at the runtime of verification, i.e., deciding if recompute needs to be triggered or not. We show how in the standard algorithmic model of computation (say, word RAM model), using suitable data structures, this can be done efficiently. Furthermore, as standard in convex optimization, we can use Nemirovski-Yudin's model of computation [5]} which counts the number of accesses to the first-order (gradient) information of the function, and a projection oracle.
Let \(\mathfrak {G}\)  denote the compute cost for one gradient access or projection in the standard model of computation –
we assume that both oracles require the same compute. In the rest of the discussion, we provide runtime as a function of the problem parameters ignoring all constants. Furthermore, since we assumed that the number of samples at any point in the stream is between \(\frac{n}{2}\)  and \(2n\) , we will just work with \(n\)  samples, and everything would still be the same, up to constants.
Verification runtime of Algorithm REF .
For Algorithm REF ,
note that for deletion, for every iteration, we need to check if the used mini-batch \(b_t\)  contained the requested point. A brute force search takes \(O(m)\)  time, whereas if we sort when we save the mini-batch indices \(b_t\) , we can do a binary search in \(O(\log {m})\)  time;
we can even do constant time search by storing a dictionary/hash table, giving us an \(O(T)\)  total time.
The most efficient way however is to store a dictionary of sample to mini-batch iterations that the sample was used in. For this, it takes \(O(1)\)  time lookup for every edit request.
For insertion, similarly, at every iteration, we first sample from a Bernoulli with bias \(m/n\)  which takes constant time, giving us \(O(T)\)  total time. However, equivalently, we just sample one Bernoulli with bias \(Tm/n\)  and recompute based on its outcome. This gives us an \(O(1)\)  time lookup for every edit request.

Verification runtime of Algorithm REF .
For Algorithm REF , we can similarly search in constant time whether the deleted point was used in any iteration or not.
For every iteration in which the deleted point is in the mini-batch, we need to compute a gradient at a new point, so as to replace the deleted point. Sampling a point uniformly from a discrete universe takes linear time (in the size) in the worst case, but with some pre-processing can be done in logarithmic/constant time. For example, when saving the mini-batch indices \(b_t\) , if we save a sorted list of the indices not sampled, using binary search, we can sample in \(O(\log {n-(m-1)})\)  time. The more efficient way is, if we save a probability table, then we can use Alias method to sample in \(O(1)\)  time [6]}. Hence for such iterations, we query two gradients, and it takes \(O(d)\)  compute to add/subtract this gradients.
Since the total number of iterations in which a deleted point was sampled in, in expectation, is \(\frac{Tm}{n}\) , the expected total compute is \(\frac{Tm(\mathfrak {G}+d)}{n}\) .
We now consider the computational cost of rejection sampling.
In Algorithm REF ,
at every iteration we check if Unif\({0,1} \le \frac{\phi _{(_{t}^{\prime },\sigma ^2)}(\xi _t)}{\phi _{(_{t},\sigma ^2)}(\xi _t)}\) , where \(\phi _{(_{t},\sigma ^2)}(\cdot )\)  and \(\phi _{(_{t}^{\prime },\sigma ^2)}(\cdot )\)  are probability densities evaluated at the sampled point \(\xi _t\) . We thus need to compute this ratio of probability densities – since these are Gaussian densities, the ratio is just the following the expression:
\(\frac{\phi _{(_{t}^{\prime },\sigma ^2)}(\xi _t)}{\phi _{(_{t},\sigma ^2)}(\xi _t)} = \frac{\frac{1}{(\sqrt{2\pi \sigma ^2})^d}\exp {-\frac{{_{t}^{\prime }-\xi _t}^2}{2\sigma ^2}}}{\frac{1}{(\sqrt{2\pi \sigma ^2})^d}\exp {-\frac{{_{t}-\xi _t}^2}{2\sigma ^2}}} = \exp {\frac{1}{2\sigma ^2}{{_{t}-\xi _t}^2 - {_{t}^{\prime }-\xi _t}^2}}.\) 
It takes \(O(d)\)  time to do the above computation. Moreover, we only need to compute the ratio in iterations where the means differ –
these correspond to the iterations where the deleted point was sampled or the inserted point would have been sampled. By a direct computation, the expected number of such iterations is \(\frac{Tm}{n}\) .
This gives us a computational cost of \(\frac{Tmd}{n}\)  for rejection sampling, and hence the expected runtime of verification is \(\frac{Tm(\mathfrak {G}+d)}{n}\) .
We now state bounds on runtime for both unlearning algorithms.
The expected total unlearning runtime of Algorithm REF  for \(k\)  edit requests is
\(O{\max {k,\min {\rho ,1} k \cdot \text{Training time}}}\) .

[Proof of claim:runtimesubsampleGD]
The total runtime of Algorithm REF  is the time for verification plus the runtime for recomputation, whenever a recompute is triggered.
The recomputation time is just the training time, and in the model considered, excepted cost of one recomputation takes \(O{T m {\mathfrak {G}+d}}\)  time, since at every iteration, \(m\)  gradients are computed and vectors added.
As discussed in sec:runtime, the expected verification time for Algorithm REF  is \(O(1)\) .
From prop:unlearn-subsample-GD, the unlearning alg:unlearn-subsample-gd recomputes with probability \(O(\min {\rho ,1})\)  for one edit request.
Therefore, using Proposition REF  which bounds the number of recomputes, we have that the expected total runtime is bounded as \(kO(1) + 4k\min {\rho ,1} \cdot O{T m {\mathfrak {G}+d}} \le O{\max {1, \min {\rho ,1}T m {\mathfrak {G}+d}}k}\) .
For a sufficiently large \(\rho \) , the unlearning time of Algorithm REF  is clearly dominated by the training time.
In particular, in the corresponding batch Algorithm REF , we set \(m=\frac{\rho n}{T}\) , giving a total runtime of \(O{\max {1, \min {\rho ^2,1}n {\mathfrak {G}+d}}k}\) .
Hence for \(\rho \gtrsim \frac{1}{\sqrt{n(\mathfrak {G}+d})}\) , the total runtime in expectation is at most \(O(\min {\rho ,1} \cdot k\cdot \text{Training time})\) .
In the other case, the expected total runtime is just \(O(k)\) .
The expected total unlearning runtime of Algorithm REF  for \(k\)  edit requests is
\(O{\max {k,\min { \rho \sqrt{T},1} \cdot k \cdot \text{Training time}}}\) .
[Proof of claim:runtimenoisyASGD]
As before, the total runtime of Algorithm REF  is the time for verification plus the runtime for recomputation, whenever a recompute is triggered.
As discussed in sec:runtime, the expected verification time for Algorithm REF  is \(O{\frac{Tm(\mathfrak {G}+d)}{n}}\) .
The recomputation in this case may be partial but it also includes a reflection.
The reflection operation with \(d\)  dimensional vectors takes \(O(d)\)  compute.
Furthermore, we upper bound the partial recomputation time by worst-case full recomputation time, giving a recomputation time
\(O{T m {\mathfrak {G}+d}}\) .
From lem:noisy-m-sgd-accept, we have that the unlearning coupling is not maximal but recomputes with probability \(\min {\rho \sqrt{T},1}\) .
Finally, by Proposition REF  we have that the expected total runtime is bounded as \(kO{T m {\mathfrak {G}+d}} + 4k\min {\rho \sqrt{T},1} \cdot O{ \frac{Tm(\mathfrak {G}+d)}{n}} \le O{\max {\min {\rho \sqrt{T},1},\frac{1}{n}}kmT(\mathfrak {G}+d)}\) .
In contrast, for Algorithm REF , the runtime is at most \(O{\max {\min {\rho \sqrt{T},1},\frac{1}{n}}kmT(\mathfrak {G}+d)}\) .
Our lower bounds will show that we need \(\rho \gtrsim \frac{1}{n}\)  to get any non-trivial accuracy.
Therefore the maximum is always obtained by \(\min {\rho \sqrt{T},1}\) .
Moreover, \(k\)  is a trivial lower bound on runtime, since we need to observe all \(k\)  edit requests.
Hence, we get that the total runtime in expectation, is at most \(O{\max {k,\min {\rho \sqrt{T},1}\cdot k \cdot \text{Training time}}}\) .

Space complexity
In this work, the objective was not to optimize the memory used, but rather, to study if the problem can be solved computationally efficiently, no matter how much (reasonable) memory the algorithm uses. However, we discuss, in this section, that the space complexities of the proposed algorithms, which we will see, is arguably, reasonably small.
We ignore the space used to store the dataset.
In both algorithms, we save a hash-table of iterations to samples - since we do \(T\)  iterations with \(m\)  samples each, this takes space of \(O(Tm)\)  words.
We also store all the iterates, which are \(d\) -dimensional vectors, so this takes a space of \(O(dT)\)  words.
Finally, we also store a dictionary of iterations to models, which takes \(O(T)\)  space. The space complexity therefore is \(O(T(\max {m,d})\) . Plugging in the values of \(T\) , we get the following.
Algorithm REF :
Plugging \(T \le \frac{\rho n}{m}\)  from prop:upperboundsubsample, we get space complexity = \(O{\rho n \max {1,\frac{d}{m}}}\) .
As remarked in sec:alg-learning-sub-sample-GD, we can improve the space complexity by not requiring to save all the iterates and yet have the same unlearning runtime. In the proof of claim:runtimesubsampleGD, we upper bound the recomputation time by a full re-computation time - this means that the upper bound on unlearning runtime holds even if the algorithm does full retraining everytime verification fails. The unlearning alg:unlearn-subsample-gd can thus be modified as follows: for deletion, instead of continue retraining from iteration \(t\)  where the deleted point participates, we can just do full retraining, with fresh randomness for all mini-batches. For insertion, note that when if condition is met (line 6 in alg:unlearn-subsample-gd), we use the iterate \(_t\)  to compute the gradient on the inserted point (line 8 in alg:unlearn-subsample-gd); however, if we don't save \(_t\) , we can just compute it on the fly by doing a full retraining with the same old mini-batches. After \(_t\)  is computed, we just continue as in alg:unlearn-subsample-gd.
With the above modification, we only need to save a hash-table of used samples to binary values which correspond to whether they were used or not, which takes \(O(Tm)\)  words, and a \(d\)  dimensional model. Hence, the space complexity of alg:unlearn-subsample-gd is \(O(Tm+d)\)  words.

Algorithm REF :
From Proposition REF , note that if \(m\le O(T^3)\) , \(T = \frac{\rho n}{md}\) , and therefore, \(dT = \frac{\rho n}{m}\) .
If we use the largest mini-batch size \(m= O(T^3)\) , then \(T = \sqrt{\frac{\rho n}{\sqrt{d}}}\) , and hence \(dT = d^{3/4}\sqrt{\rho n}\) .
Therefore, the space complexity is \(O(T\max {m,d}) \le O{\max {\frac{(\rho n)^2}{d},d^{3/4}\sqrt{\rho n}}}\)  words.

Other algorithms and batch unlearning
To demonstrate the generality of our framework, we give two more algorithms. The first is noisy-m-SGD which is the same as alg:noisy-m-a-sgd but without acceleration, and the second is quantized-m-SGD, based on randomized quantization. We note that both algorithms have worse theoretical guarantees than alg:noisy-m-a-sgd, however the first establishes our claim that acceleration is beneficial, whereas the second shows how a previous work of [7]} for \(k\) -means clustering, can, not only be seen as a special case of our framework, but also extended to general convex risk minimization problems.
Moreover, in the second case, we consider a more general setup of batch edit requests, and show that our techniques are flexible enough to easily generalize to the batch variant.
noisy-m-SGD
[ht]
noisy-m-SGD(\(_{t_0},t_0\) )

[1]
Initial model \(_{t_0}\) , data points \({_1,\ldots , _n}\) , \(T, \eta , m\) 
\(_0 = 0\) 
\(t=t_0,t_0+1,\ldots , T\) 
Sample mini-batch \(b_t\)  of size \(m\)  uniformly randomly
\(_t = \frac{1}{m}\sum _{j \in b_t} \nabla f(_t,_j) \) 
Sample \(\theta _t \sim (0,\sigma ^2_d)\) 
\(_{t+1} = {_t - \eta { _t +\theta _t}}\) 
Save(\(b_t, \theta _t,_t, _t\) )
\(\hat{}_S = \frac{1}{T}\sum _{t=1}^{T+1}_t\)
Let \(f(., )\)  be an \(L\) -smooth \(G\) -Lipschitz convex function \(\forall \ \) .
Algorithm REF , run with \(t_0=1, \eta = \min {\frac{1}{2L}, \frac{D}{{\frac{G}{\sqrt{m}} +\sigma }\sqrt{T}}}\) , \(\sigma = \frac{8\sqrt{T}G}{n\rho }\) , and \(T\ge \frac{(\rho n)^2}{16m^2}\)  outputs \(\hat{}_S\)  which is \(\min {1,\rho }\) -TV stable and satisfies \({\hat{F}_S(\hat{}_S) - \hat{F}_S(^*_S)} \lesssim \frac{GD\sqrt{d}}{\rho n}\) 

[Proof of prop:upperboundnoisy-m-sgd]
The \(TV\) -stability guarantee of \(\rho = \frac{8\sqrt{T}G}{n}\)  follows exactly as in the proof of prop:upper-bound-noisy-m-a-sgd.
We now proceed to the accuracy guarantee, which follows simply by guarantee of SGD on smooth convex functions. We have already shown in prop:upper-bound-noisy-m-a-sgd that the gradients are unbiased and its variance bounded by \(\frac{2G^2}{m} + \frac{2G^2}{m} + \sigma ^2d\) .
Therefore, using Theorem 4.1 in [8]} with step-size \(\eta \le \frac{1}{2L}\) , we have
\({ \hat{F}(\hat{}) - \hat{F}(^*)} \le O{ 2\eta ^2 + \frac{D^2}{\eta T}} = O{2\eta {\frac{2G^2}{m} +\sigma ^2d} + \frac{D^2}{\eta T}}\) 
Let \(\tilde{G}^2 = \frac{2G^2}{m} +\sigma ^2d\) , balancing the trade-off in \(\eta \)  gives us \(\eta = \frac{D}{\tilde{G}\sqrt{T}}\) . Therefore setting \(\eta = \min {\frac{1}{2L},\frac{D}{\tilde{G}\sqrt{T}}}\)  gives us
\({ \hat{F}(\hat{}) - \hat{F}(^*)}&\le O{\frac{LD^2}{T} + \frac{\tilde{G}D}{\sqrt{T}}} \le O{\frac{LD^2}{T} + \frac{GD}{\sqrt{Tm}} + \frac{\sigma \sqrt{d} D}{\sqrt{T}}} \\&\le O{\frac{LD^2}{T} + \frac{GD}{\sqrt{Tm}} + \frac{G D\sqrt{d}}{n \rho }}\) 
Finally, the condition in the sub-sampling amplification \(\frac{8G^2}{m^2\sigma ^2} \le 1.256\)  again becomes \(T\ge \frac{(n\rho )^2}{16m^2}\) .
We now show that alg:noisy-m-sgd achieves the same upper bound on excess empirical risk as alg:noisy-m-a-sgd.
Let \(f(., )\)  be an \(L\) -smooth \(G\) -Lipschitz convex function \(\forall \ \) .
Algorithm REF , run with \(m \ge \min {\frac{d}{16}, \frac{1}{4}{\frac{(\rho n) G\sqrt{d}}{LD}}^{1/2}}\) , \(\eta = \min {\frac{1}{2L}, \frac{D}{{\frac{G}{\sqrt{m}} +\sigma }\sqrt{T}}}\) , \(\sigma = \frac{8\sqrt{T}G}{n\rho }\) , and \(T=\max {\frac{(\rho n)^2}{md}, \frac{LD\rho n}{G\sqrt{d}}}\)  outputs \(\hat{}_S\)  which is \(\rho \) -TV stable and satisfies \({\hat{F}_S(\hat{}_S) - \hat{F}_S(^*_S)} \lesssim \frac{GD\sqrt{d}}{\rho n}\)
[Proof of cor:upper-bound-noisy-m-sgd]
We start with the result in prop:upperboundnoisy-m-sgd.
Note that as long as \(\frac{GD}{\sqrt{mT}} \gtrsim \frac{LD^2}{T} \iff m \lesssim \frac{TG^2}{LD}\) , the second term is larger than the first.
We balance the two trade-offs in \(T\) .
Optimizing the trade-off between second and third term gives us
\(\frac{GD}{\sqrt{mT}} = \frac{GD\sqrt{d}}{\rho n} \iff T = \frac{(\rho n)^2}{md}\) ;
and optimizing the second trade-off gives us
\(\frac{\sqrt{d}}{\rho n} = \frac{LD^2}{T} \iff T = \frac{LD^2(\rho n)}{\sqrt{d}}\) . Hence setting \(T = \max {\frac{(\rho n)^2}{md},\frac{LD^2(\rho n)}{\sqrt{d}} }\)  yields an expected excess empirical risk of \(O{\frac{GD\sqrt{d}}{n \rho }}\) .
We now look at the condition \(T\ge \frac{(n\rho )^2}{16m^2}\)  given in prop:upperboundnoisy-m-sgd, with \(T\)  set as \(T = \max {\frac{(\rho n)^2}{md},\frac{LD^2(\rho n)}{\sqrt{d}} }\) .
We therefore require \(\frac{(\rho n)^2}{md} \ge \frac{(\rho n)^2}{16m^2} \iff m \ge \frac{d}{16}\) , as well as \(\frac{LD(\rho n)}{G\sqrt{d}} \ge \frac{(\rho n)^2}{16m^2} \iff m \ge \frac{1}{4}{\frac{(\rho n) G\sqrt{d}}{LD}}^{1/2}\)  - this recovers the condition \(m \ge \min {\frac{d}{16}, \frac{1}{4}{\frac{(\rho n) G\sqrt{d}}{LD}}^{1/2}}\)  in the Proposition statement.
Hence, combining all the above arguments, we get that for any \(m \ge \min {\frac{d}{16}, \frac{1}{4}{\frac{(\rho n) G\sqrt{d}}{LD}}^{1/2}}\)  , setting \(T = \max {\frac{(\rho n)^2}{md},\frac{LD(\rho n)}{G\sqrt{d}}}\) , yields an expected excess empirical risk of \(O{\frac{GD\sqrt{d}}{n \rho }}\) .
[ht]
Unlearning for noisy-m-SGD

[1]
Delete point with index \(j\)  or insert \(\)  (with index \(n+1\) ) for noisy-m-SGD
\(t=1,2\ldots , T\) 
Load\({\theta _t, _t, b_t, _t)}\) 
deletion and \(j\in b_t\) 
Sample \(i \sim \text{Uniform}([n]\backslash b_t)\) 
\( _t^{\prime } = _t - \frac{1}{m}{\nabla f(_t,_j) - \nabla f(_t,_i)}\) 
Save(\( _t^{\prime }, b_t\backslash {j} \cup {i} \) )
insertion and Bernoulli\({\frac{m}{n+1}}\) 
Sample \(i \sim \text{Uniform}(b_t)\) 
\(_t^{\prime } = _t- \frac{1}{m}{\nabla f(_t,_i) - \nabla f(_t,)}\) 
Save(\( _t^{\prime }, b_t\backslash {i} \cup {n+1} \) )
\(\xi _t = _t + \theta _t\) 
Uniform\({0, 1} \ge \frac{\phi _{(_t^{\prime },\sigma ^2)}(\xi _t)}{ \phi _{(_t,\sigma ^2)}(\xi _t)}\) 
\( \xi ^{\prime }_t = \text{reflect}(\xi _t,_t^{\prime }, _{t})\) 
\(_{t+1} = _t - \eta \xi _t^{\prime }\) 
\(\text{Save}(\xi ^{\prime }_t)\) 
noisy-m-SGD(\(_{t+1},t+1\) )  // Continue retraining, on current dataset
The choice of \(T\)  in Proposition REF  yields that the largest mini-batch size that can be set, without hurting runtime, is \(m = \frac{\rho n G}{\sqrt{d} L D}\) .
Furthermore, the condition \(m \ge \min {\frac{d}{16}, \frac{1}{4}{\frac{(\rho n) G\sqrt{d}}{LD}}^{1/2}}\)  becomes \(T\ge {\frac{\sqrt{d}LD}{4G}}^2\) .
We now state and prove the main theorem for this section.


Let \(f(\cdot , )\)  be an \(L\) -smooth \(G\) -Lipschitz convex function \(\forall \ \) .
For any \(\frac{1}{n} \le \rho \le 1\) , using alg:noisy-m-sgd as the learning algorithm and alg:unlearn-noisy-m-sgd as its unlearning algorithms, then given a stream of edit requests,

Satisfies exact unlearning at every point in the stream.

At time \(i\)  in the stream of edit requests, outputs \(\hat{}_{S^i}\) , such that if \({\frac{L^{1/2}D^2 \sqrt{d}}{G (\rho n)}}^{2/3} \le \frac{GD}{\sqrt{\rho n}}\) , then its with excess empirical risk bounded as,
\({\hat{F}_S(\hat{}_{S^i}) - \hat{F}_{S}(_{S^i}^*)} \lesssim {\frac{L^{1/2}D^2 \sqrt{d}}{G (\rho n)}}^{2/3}\) 

For \(k\)  edit requests, the expected total unlearning runtime is \(O(\max {\rho k \cdot \text{Training time},k})\)

[Proof of thm:main-thm-noisy-sgd]
We proceed as in the proof of thm:main-result. For any \(0< \tilde{\rho }\le 1\) , from prop:upperboundnoisy-m-sgd, the output \(\hat{}_{S}\)  is \(\tilde{\rho }\) -TV stable, and the excess empirical risk using alg:noisy-m-sgd on a dataset \(S\)  on \(n\)  points, is bounded as,
\(\mathbb {E}\hat{F}_S(\hat{}_S) - \hat{F}_S(^*_S) \le O{\frac{LD^2}{T} + \frac{GD}{\sqrt{Tm}} + \frac{G D\sqrt{d}}{n \tilde{\rho }}}\) 
It can be easily verified that prop:unlearn-noisy-m-A-SGD and prop:unlearn-noisy-m-A-SGD still holds for noisy-m-SGD, which gives us that the algorithm satisfies exact unlearning at every time in the stream, proving the first part of the claim, Moreover, its recompute probability bounded by \(O(\tilde{\rho }k \sqrt{T})\)  and therefore the unlearning runtime bounded by \(O(\max {k,\tilde{\rho }k \sqrt{T}\cdot \text{Training time}}\) . Substituting \(\tilde{\rho }= \frac{\rho }{\sqrt{T}}\) , and using the largest mini-batch size \(m = {\frac{G}{LD}}^2T\) , the upper bound on excess empirical risk becomes \(\frac{LD^2}{T} + \frac{G D\sqrt{d}\sqrt{T}}{n \rho } \) . Optimizing the trade-off, we have \(\frac{LD^2}{T} = \frac{LD^2}{T} \iff T = {\frac{LD (\rho n)}{G\sqrt{d}}}^{2/3}\) , and the excess empirical risk bound upper bound is \(\frac{LD}{T} = {\frac{L^{1/2}D^2 \sqrt{d}}{G (\rho n)}}^{2/3}\) .
Note that this also proves the third part of the claim. Furthermore, as in the proof of thm:main-result, it can be verified that the condition \(T\ge \frac{(\tilde{\rho }n)^2}{16m^2}\)  is equivalent to \({\frac{L^{1/2}D^2 \sqrt{d}}{G (\rho n)}}^{2/3} \le \frac{1}{\sqrt{\rho n}}\) , which just means that the excess empirical risk of noisy-m-SGD is at most that of sub-sample-GD.
Finally, the upper bound holds for any point in the stream using the assumption that the number of samples are between \(\frac{n}{2}\)  and \(2n\) , thereby establishing the second claim.

quantized-m-SGD
The work of [7]} considers unlearning in \(k\) -means clustering. The key algorithmic technique is randomized quantization of vectors to a \(\tau \) -lattice.
The intuition is that if the vector is an average of \(n\)  data points which are bounded in norm, then upon changing one data point, the vectors \(O{\frac{1}{n}}\)  close.
Therefore, if the lattice is sufficiently coarse, then it would ensure that both are mapped up the same point in the lattice. However, if we consider deterministic quantization, then there exists points such that for any \(\epsilon >0\) , shifting the point by \(\epsilon \)  changes the quantized point. Therefore, we first shift the lattice by a uniformly random phase, which ensures that such a situation occurs with a small probability.
In their application of \(k\) -means clustering, this vector is a cluster centroid, which is an average of the data points in the cluster. We apply this idea to convex risk minimization problems, wherein we quantize average gradients, which by the Lipshcitzness assumption are bounded in norm.
We now introduce the quantization operation formally.
Given a vector \(\) , let \(\theta \sim \text{Unif}\left[-\frac{1}{2},-\frac{1}{2}\right]^d\) , consider the quantization given by:
\(Q_{\theta }() = \tau {\theta + \arg \min _{j \in ^d} {- \tau (\theta + j)}}\) 
We now state a result about the quantization operation.
Let \(B_\delta (\)  denote the Euclidean call of radius \(\delta \)  centered at \(.The following holds for the quantization operation,\begin{enumerate}\item For any , {Q_\theta ()}=  and {{-{}}}^2 \le \tau ^2d \item For any vector , ¶{ \exists B_\delta (:Q_\theta (\ne _\theta (} \le \frac{2d \delta }{\tau }\end{enumerate}\)
[Proof of claim:quantize]
Note that for a given \(\) , \(_\theta () \sim \text{Unif}\left[-\frac{\tau }{2},+\frac{\tau }{2}\right]^d\) , hence \({Q_\theta ()}= \) . Furthermore, since \(-{} \sim \text{Unif}\left[-\frac{\tau }{2},\frac{\tau }{2}\right]^d\) , we have \({{-{}}}^2 =d \mathbb {E}{_1-{_1}}^2 = \frac{d\tau ^2}{12}\) .
The second part of the claim is Lemma C.2 in [7]}.
To see why [7]} is a special case of our framework, note that the total variation distance between two random variables is at most the probability of disagreement under any coupling. [7]} uses the same quantization randomness (used for training) for verifying after the edit request - this corresponds to a trivial coupling between the quantization randomness, hence the total variation distance between the outputs is bounded by the upper bound on the probability that the quantized points change (see claim:quantize). This establishes that it is a TV stable method. Finally, as said before, using the same quantization randomness corresponds to a trivial coupling, but can be shown to be maximal since the probability distribution is uniform around the to-be-quantized point. Therefore, we have that [7]} uses a maximal coupling based unlearning method.
Batch unlearning:
We consider a batch unlearning setup, wherein instead of observing an insertion or deletion request, we observe a batch edit request with insertions and deletions. We demonstrate that our general approach of coupling mini-batch indices is flexible enough to handle this variant naturally. The batch unlearning ideas and results extend to other algorithms: noisy-m-A-SGD, noisy-m-SGD and subsample-GD. We also note that the computational benefit of batch unlearning as opposed to handling edits one by one is only a constant factor, which at best is two.
We now discuss how we extend the randomized quantization idea to convex risk minimization.
In our learning algorithm quantized-m-SGD, at each iteration, we draw a mini-batch of \(m\)  samples, uniformly randomly from \(n\)  samples, use it to compute the gradient on the previous iterate , quantize using a randomly sampled phase, and update. Algorithm REF  implements the above procedure.
[ht]
quantized-m-SGD

[1]
Initial model \(_1\) , data points \({_1,\ldots , _n}\) ,\(T,\eta \) 
\(t=1,2\ldots , T\) 
Sample mini-batch \(b_t\)  of size \(m\)  uniformly randomly
Sample \(\theta _t \sim \text{Unif}\left[-\frac{1}{2},\frac{1}{2}\right]^d\) 
\(_t = \frac{1}{m}\sum _{j \in b_t} \nabla f(_t,_j)\) 
\(_{t+1} = _t - \eta Q_{\theta _t}{_t}\) 
Save(\(\theta _t,_t, b_t, _t\) )
\(\hat{}_S = \frac{1}{T+1}\sum _{j=1}^{T+1} _j\)
We first prove a lemma which bounds the total variation distance between outputs generated by quantized-m-SGD on arbitrarily differing datasets - these can be thought of as arising after a batch edit request.
Let \(S\)  and \(S^{\prime }\)  be two datasets of \(n\)  and \(n+k_2\)  points respectively, such that \(S\)  has \(k_1\)  points which differ from \(S^{\prime }\)  i.e. \({S\backslash S^{\prime }}=k_1\) , therefore \(S\)  and \(S^{\prime }\)  differ by \(k_1+k_2\)  points.
Let \({_j}_{j=1}^T\)  and \({_j^{\prime }}_{j=1}^T\)  be iterates of quantized-m-SGD on datasets \(S\)  and \(S^{\prime }\)  respectively. The total variation distance between distribution of average iterates \(\hat{}_S\)  and \(\hat{}_{S^{\prime }}\)  is bounded as,
\(\text{TV}(\hat{}_S, \hat{}_{S^{\prime }}) \le \frac{4GTd(k_1+k_2)}{n\tau }\) 
[Proof of lemma:tv-stability-qsgd]
Without loss of generality, we enumerate \(S\)  and \(S^{\prime }\)  into subsets as follows: let \(S_1\)  and \(S_1^{\prime }\)  be the first \(n-k_1\)  elements of \(S\)  and \(S^{\prime }\)  which are the same. Let \(S_2\)  and \(S_2^{\prime }\)  be the next \(k_1\)  differing elements in \(S\)  and \(S^{\prime }\)  respectively. Finally, let \(S_3^{\prime }\)  be the last \(k_2\)  elements of \(S^{\prime }\) .
We look at iteration \(t\)  of quantized-m-SGD and fix the previous model \(_t = _t^{\prime } = \) . We will now compute the conditional total variation distance between \(_{t+1}\)  and \(_{t+1}^{\prime }\) . Note that since the only randomness is in the sub-sampling and quantization, we can compute the total variation distance between sub-sampled quantized gradients on fixed \(\)  for both datasets, and this will lower bound total variation distance between the iterates \(_{t+1}\)  and \(_{t+1}^{\prime }\)  by data processing inequality.
Let \(b^{(1)}\)  and \(b^{(2)}\)  be a uniform sample of \(m\)  points from datasets \(S\)  and \(S^{\prime }\)  respectively.
For a fixed \(\) , let the gradient on \(S\)  indexed by \(b^{(1)}\)  be denoted as \(_{b^{(1)}}^S() = \frac{1}{m}\sum _{j \in b^{(1)}} \nabla f(,_{j}^S)\) , and similarly for \(S^{\prime }\) .
Let \(P\)  and \(Q\)  denote the probability distribution of \(_{b^{(1)}}^S()\)  and \(_{b^{(2)}}^{S^{\prime }}()\)  respectively.
We have the following claim, which we will prove via mathematical induction on \(k_2\) : for any measurable set \(R\) , for any \(k_2\) , \( {P(R)-Q(R)} \le \frac{4Gdk_1}{n\tau } + \sum _{j=1}^{k_2}\frac{4Gd}{(n+j)\tau }\) .

Base case 1: \(\mathbf {k_2=0:}\)
Firstly note that both \(b^{(1)} \sim \text{Unif}([n],m)\)  and \(b^{(2)} \sim \text{Unif}([n],m)\) , and consider the trivial coupling \(b^{(1)}=b^{(2)}=b\) , where \(b \sim \text{Unif}([n],m)\) , be a uniform sample of \(m\)  points from \([n]\) .
We now use the fact that total variation distance is at most the probability of disagreement for any coupling. This gives us that
\(\text{TV}(Q_\theta (_{b^{(1)}}^S()), Q_\theta (_{b^{(2)}}^{S^{\prime }}()) \le ¶{Q_\theta (_b^S()) \ne Q_\theta (_b^{S^{\prime }}())}\) 

We will focus on upper bounding the right hand side.
The proof follows by using the quantization guarantee (Claim REF ) combined amplification from subsampling.
Without loss of generality, assume that the first \(k_1\)  samples in \(S\)  and \(S^{\prime }\)  are the ones that differ. Fix the random (uniform) sample \(b\)  of indices - suppose for this fixed value of \(b\) , exactly \(j\)  differing data points are sampled. From \(G\)  Lipschitzness, and that we have exactly \(j\)  differing data points, \({_b^S() - _b^{S^{\prime }}() } \le \frac{2Gj}{m}\) . Hence, applying Claim REF , we have that
\(¶{Q_\theta (_b^S()) \ne Q_\theta (_b^{S^{\prime }}()) \Big \vert b \text{ producing } j \text{ differing samples}} \le \frac{4 Gdj}{\tau m}.\) 
We will now integrate with respect to the randomness in \(b\)  - for this, we need to calculate the probability that a sample of \(b\)  (uniform \(m\)  out of \(n\) ) produces exactly \(j\)  differing data points, call it \(p(j)\) . By direct computation, we have that \(p(j) = \frac{{k_1 \atopwithdelims ()j}{n-k_1 \atopwithdelims ()m-j}}{{n \atopwithdelims ()m}}\) . Hence we have,
\(¶{Q_\theta (_b^S()) \ne Q_\theta (_b^{S^{\prime }}())} &= \sum _{j=0}^{k_1} p(j) ¶{Q_\theta (_b^S()) \ne Q_\theta (_b^{S^{\prime }}()) \Big \vert \ b \ \text{produces } j \ \text{differing samples}} \\& \le \sum _{j=0}^{k_1} \frac{{{k_1} \atopwithdelims ()j}{n-{k_1} \atopwithdelims ()m-j}}{{n \atopwithdelims ()m}} \frac{4 Gdj}{\tau m}= \frac{m{k_1}}{n} \frac{4 Gd}{\tau m} = \frac{4 Gd{k_1}}{\tau n}\) 

where the second last equality is a consequence of Vandermonde's identity, as we show below. We need to show that \(\sum _{j=0}^{k_1} \frac{{{k_1} \atopwithdelims ()j}{n-{k_1} \atopwithdelims ()m-j}j}{{n \atopwithdelims ()m}} = \frac{m{k_1}}{n} \iff \sum _{j=0}^{k_1} {{k_1} \atopwithdelims ()j}{n-{k_1} \atopwithdelims ()m-j}j = \frac{m{k_1}}{n} {n \atopwithdelims ()m} = {k_1}{n-1 \atopwithdelims ()m-1}\) . This holds because,
\(\sum _{j=0}^{k_1} {{k_1} \atopwithdelims ()j}{n-{k_1} \atopwithdelims ()m-j}j &= \sum _{j=0}^{k_1} \frac{{k_1}}{j}{{k_1}-1 \atopwithdelims ()j-1}{n-{k_1} \atopwithdelims ()m-j}j = {k_1}\sum _{j=0}^{k_1}{{k_1}-1 \atopwithdelims ()j-1}{n-{k_1} \atopwithdelims ()m-j} \\&={k_1}\sum _{j=0}^{{k_1}-1}{{k_1}-1 \atopwithdelims ()j}{n-{k_1} \atopwithdelims ()(m-1)-j} = {k_1}{n-1 \atopwithdelims ()m-1}\) 

where in the second last equality, we re-indexed the sum which removes the first element, but it was zero anyway, and the last equality follows from Vandermonde's identity.

Base case 2: \(\mathbf {k_2=1:}\)
In this case, \(S^{\prime }\)  has one more element that \(S^{\prime }\)  - let this point be denoted as \(q\) .
In this case, the probability distribution using \(S^{\prime }\)  has the form \(Q = {1-\frac{m}{n+1}}Q_1+\frac{m}{n+1}Q_2\) , where \(Q_1\)  is the probability distribution conditioned on the event that \(q\)  is sub-sampled, and \(Q_2\)  is the probability distribution conditioned on the complementary event. For any measurable set \(R\) , we have,
\({P(R)-Q(R)} = {{1-\frac{m}{n+1}}Q_1(R)+\frac{m}{n+1}Q_2(R)-P(R)}\) 

Note that \(Q_1,Q_2\)  and \(P\)  are all probability distributions over \(n\)  elements. Furthermore, \(P\)  and \(Q_1\)  are probability distributions over \(k_1\)  differing elements, therefore we can use base case \(k_2=0\)  to get that \({P(R)-Q_1(R)} \le \epsilon _1:=\frac{4Gdk_1}{n\tau }\) . We therefore get,
\({P(R)-Q(R)} &\le {{1-\frac{m}{n+1}}Q_1(R)+\frac{m}{n+1}Q_2(R)-Q_1(R)+\epsilon _1} \\& = {\frac{m}{n+1}{Q_2(R)-Q_1(R)} + \epsilon _1}\) 

Finally note that \(Q_1\)  and \(Q_2\)  are probability distributions over \(n\)  such that upon sub-sampling \(m\)  elements, there is exactly one differing element, therefore
we get, \({Q_2(R)-Q_1(R)} \le \frac{4Gd}{m\tau }\) . We therefore have that
\({P(R)-Q(R)} &\le \frac{m}{n+1}\frac{4Gd}{m\tau } + \frac{4Gdk_1}{n\tau } = \frac{4Gd}{(n+1)\tau } + \frac{4Gdk_1}{n\tau }\) 

Induction Hypothesis:
Suppose the following holds for \(k_2 \le \tilde{k}\) : for any measurable set \(R\) , \({P(R)-Q(R)} \le \frac{4Gdk_1}{n\tau } + \sum _{j=1}^{\tilde{k}}\frac{4Gd}{(n+j)\tau }\) .

Induction Step: \(k_2 = \tilde{k}+1\)
Let the last element of \(S^{\prime }\)  be \(q\) . As in the base case, we decompose the distribution \(Q\)  into a mixture of two components based on whether \(q\)  is sampled or not. We have \(Q = {1- \frac{m}{n+\tilde{k}+1}} Q_{1} + \frac{m}{n+\tilde{k}+1}Q_2\) . Note that \(Q_1\)  is a probability distribution which does not use the last element of \(S^{\prime }\) . Therefore we can use Induction hypothesis which gives us that \({Q_1(R)-P(R)}\le \frac{4Gdk_1}{n\tau } + \sum _{j=1}^{\tilde{k}}\frac{4Gd}{(n+j)\tau }\) . We therefore get,
\({P(R)-Q(R)} &= {P(R)-{1- \frac{m}{n+\tilde{k}+1}} Q_{1}(R) - \frac{m}{n+\tilde{k}+1}Q_2(R)} \\& \le {\frac{m}{n+\tilde{k}+1}{Q_1(R)-Q_2(R)} + \frac{4Gdk_1}{n\tau } + \sum _{j=1}^{\tilde{k}}\frac{4Gd}{(n+j)\tau }} \\&\le {\frac{m}{n+\tilde{k}+1} \frac{4Gd}{m\tau } + \frac{4Gdk_1}{n\tau } + \sum _{j=1}^{\tilde{k}}\frac{4Gd}{(n+j)\tau }} \\& = \frac{4Gdk_1}{n\tau } + \sum _{j=1}^{\tilde{k}+1}\frac{4Gd}{(n+j)\tau }\) 

where in the last inequality, as in the base case, we used that fact that distributions \(Q_1\) 
and \(Q_2\)  differ because in one we subsample the last element where as in the other we don't, so from Claim REF , for two data sets of size \(m\)  differing in one element, the failure probability is \(\frac{4Gd}{m\tau }\) . This completes the induction argument.
We bound the sum simply as \(\sum _{j=1}^{k_2}\frac{4Gd}{(n+j)\tau } \le \frac{4Gdk_2}{n\tau }\) , which gives us that the whole term is bounded by \(\frac{4Gd(k_1+k_2)}{n\tau }\) .
The above, by an application of data processing inequality, shows that the conditional TV distance between \(_{t+1}\)  and \(_{t+1}^{\prime }\)  is at most \(\frac{4Gd(k_1+k_2)}{n\tau }\) . Note that the upper bound holds uniformly over all conditioning events.
Moreover, from the maximal coupling characterization of TV distance, we have that for any coupling of \(_{t+1}\)  and \(_{t+1}^{\prime }\) ,
the conditional probability of disagreement is at most \(\frac{4Gd(k_1+k_2)}{n\tau }\) . Consider the coupling which just concatenates all these couplings, then an application of union bound over the \(T\)  iterates, the joint probability of disagreement under this coupling is at most \(\frac{4Gd(k_1+k_2)T}{n\tau }\)  which gives us our upper bound on TV distance between joint iterates.
Finally, by data processing inequality, the same upper bound holds for the average iterates which finishes the proof.
We now establish the guarantees on the learning alg:qsgd.
To handle batch edit request, we extend the notion of exact unlearning with one edit request to batch request: we term it exact batch unlearning. We similarly also extend \(\rho \) -TV-stability to \((k_1,k_2,\rho )\) -TV stability, which is \(\rho \) -TV stability under arbitrary \(k_1\)  deletions and \(k_1+k_2\)  insertions, as well as \(k_1\)  insertions and \(k_1+k_2\)  deletions.


Let \(f(., )\)  be an \(L\) -smooth \(G\) -Lipschitz convex function \(\forall \ \) .
Algorithm REF , run with \(\eta = \min {\frac{1}{2L}, \frac{D}{{\frac{G}{\sqrt{m}} +\tau \sqrt{d}}\sqrt{T}}}\) , \(\tau = \frac{4GdT}{\rho n}\) , and \(T = \max {\frac{(\rho n)}{d^{3/2}\sqrt{m}},{\frac{LD (\rho n)}{Gd^{3/2}}}^{2/3} }\)  outputs \(\hat{}_S\)  which is \((k_1,k_2,(k_1+k_2)\rho )\) -TV stable and satisfies \( { \hat{F}_S(\hat{}_S) - \hat{F}_S(^*_S)} \lesssim {\frac{\sqrt{L}GD^2 d^{3/2}}{\rho n}}^{2/3}\) .
[Proof of prop:upper-bound-qsgd]
The \((k_1,k_2,(k_1+k_2)\rho )\) -TV stability guarantee follows from lemma:tv-stability-qsgd by taking a supremeum over all datasets \(S\)  and \(S^{\prime }\)  of sizes \(n\)  and \(n+k_2\)  (or \(n-k_2)\)  to get that that TV stability is uniformly upper bound by \(\frac{4GTd(k_1+k_2)}{n\tau } =(k_1+k_2)\rho \) , where the equality follows upon setting \(\tau = \frac{4GdT}{\rho n}\) .
For the excess empirical risk bound, we use the guarantee on excess empirical risk of SGD on smooth convex functions (for example, Theorem 4.1 from [8]}), combined with the fact in claim:quantize that quantization produces unbiased estimates of the gradient with bounded variance \(^2 \le \frac{2G^2}{m} + \tau ^2 d = \frac{2G^2}{m} +\frac{16G^2d^3T^2}{(\rho n)^2}\) . Therefore, choosing step size \(\eta \le \frac{1}{2L}\) , we get
\({ \hat{F}_S(\hat{}_S) - \hat{F}_S(^*_S)} \le O{ 2\eta ^2 + \frac{D^2}{\eta T}} = O{2\eta {\frac{2G^2}{m} +\frac{16G^2d^3T^2}{(\rho n)^2}} + \frac{D^2}{\eta T}}\) 
Define \(\tilde{G}^2 = 2 {\frac{2G^2}{m} +\frac{16G^2d^3T^2}{(\rho n)^2}} \)  and set \(\eta = \min {\frac{1}{2L}, \frac{D}{\tilde{G}\sqrt{T}}}\) , which makes the upper bound
\({ \hat{F}_S(\hat{}_S) - \hat{F}(^*_S)} &\le O{ 2\eta ^2 + \frac{D^2}{\eta T}} = O{\frac{LD^2}{T} + \frac{\tilde{G}D}{\sqrt{T}}}\\&\le O{\frac{LD^2}{T}+\frac{GD}{\sqrt{Tm}} + \frac{GD d^{3/2}\sqrt{T}}{(\rho n)}}\) 
Balancing the trade-off between the las two terms gives us \(T=\frac{(\rho n)}{d^{3/2}\sqrt{m}}\) . Similarly, balancing the trade-off between the first and last term gives us \(T = {\frac{LD (\rho n)}{Gd^{3/2}}}^{2/3}\) . Hence setting \(T = \max {\frac{(\rho n)}{d^{3/2}\sqrt{m}},{\frac{LD (\rho n)}{Gd^{3/2}}}^{2/3} }\)  gives us that the expected excess empirical risk is bounded by \({\frac{\sqrt{L}GD^2 d^{3/2}}{\rho n}}^{2/3}\)  and completes the proof.
We see that the TV stability parameter above is \((k_1+k_2)\rho \)  as opposed to \((k_1+2k_2)\rho \)  which is what we would obtain with \(\rho \) -TV stability for one edit request and using the triangle inequality of TV distance (see remark:tvuptok).
The largest mini-batch size, without hurting runtime, is \(m = {\frac{G}{LD}}^2T = {\frac{G^2 (\rho n)}{(LD)^2d^{3/2}}}^{2/3}\) , which gives us \(T={\frac{LD (\rho n)}{Gd^{3/2}}}^{2/3}\) .
We now proceed to unlearning. The unlearning algorithm (alg:unlearnqsgd) upon observing an edit request comprising of both insertions and deletions, first couples the mini-batch indices (described formally in the next paragraph), and computes the gradient on the new mini-batch It then uses the same quantization randomness as in training, and checks if the quantized point changes. If it does, in any iteration, then it calls recompute.
The use of the same quantization randomness corresponds to a trivial coupling between the quantization randomness.
We explain the coupling procedure in detail below.

Batch coupling: 
We setup some notation. Consider the training dataset \(S\)  and dataset realized after the batch edit request \(S^{\prime }\) .
Given a vector \(\) , let \(Q_{\theta _1}(_{b_1}^S())\)  denote the quantized gradient vector where \(\theta _1\)  is the quantization randomness and \(b_1\)  is the mini-batching randomness on dataset \(S^{\prime }\) . Similarly, \(Q_{\theta _2}(_{b_1}^{S^{\prime }}())\)  denotes the quantized vector with \(\theta _2\)  as the quantization randomness and \(b_2\)  as the mini-batching randomness on dataset \(S^{\prime }\) . We couple \(\theta _1\)  and \(\theta _2\)  by considering the trivial coupling \(\theta _1 = \theta _2\)  i.e. the joint probability measure is defined only on the diagonal of the product measure. To couple the mini-batch indices, we consider two cases: if the training dataset \(S\)  has less more or more points than \(S^{\prime }\) . For simplicity, alg:unlearnqsgd is the pseudo-code corresponding only to the first case.
In the first case, suppose \(S\)  has \(n\)  points and \(S^{\prime }\)  has \(n+k_2\)  points, realized after \(k_1 \)  deletions and \(k_1+k_2\)  insertions. Without loss of generality, order the two datasets as follows: the first \(n-k_1\)  points in \(S\)  and \(S^{\prime }\)  are the same, call these \(S_1=S_1^{\prime }\) , next we have the last \(k_1\)  points of \(S\) , and arbitrary \(k_1\)  points of \(S^{\prime }\)  - call these \(S_2\)  and \(S^{\prime }_2\) , and moreover let the mapping of indices from \(S_2 \rightarrow S_2^{\prime }\)  by denoted by \(\iota \) . Finally we have the rest of \(k_2\)  points of \(S^{\prime }\) , call this \(S_3^{\prime }\) .
In the following discussion, and in alg:unlearnqsgd, when we consider elements of these sets, we mean their indices.
As before, let \(\mu _{n,m}\)  and \(\mu _{n+k_2,m}\)  denote the probability measures correspondingly to sampling \(m\)  elements uniformly from a discrete universe of size \(n\)  (i.e. \(S\) ) and \(n+k_2\)  (i.e. \(S^{\prime }\) ) respectively.
These sub-sampling measures are coupled in the following way in alg:unlearnqsgd.
We first sample \(b^{(1)}\sim \mu _{n,m}\)  (during training). Let \(b^{(1)}_2\)  be the \(m_2\)  indices in \(S_2\) : replace these by the corresponding indices in \(S_2^{\prime }\)  i.e. \(\iota (b^{(1)} \cap S_2)\) .
Next, sample \(b\sim \mu _{n+k_2,m}\) : let \(b_3\)  be the \(m_3\)  indices which are in \(S_3^{\prime }\) . We now resample \(b^{(2)}_1 \sim \text{Unif}(b^{(1)},m_3)\)  - these are indices used in training, which are now to be replaced. Define \(b^{(2)} = {b^{(1)} \backslash {b^{(1)} \cap S_2} \cup {\iota (b^{(1)} \cap S_2)}}\backslash b^{(2)}_1 \cup {S_3^{\prime } \cap b}\) . Let the distribution of \(b^{(2)}\)  produced in the above way be denoted as \(\mu ^{\text{edit}}_{n, m}\) .
We now show that \((b^{(1)}, b^{(2)})\)  is indeed a coupling of \(\mu _{n,m}\)  and \(\mu _{n+k_2,m}\) .
With the construction described above, we have that \(b^{(1)} \sim \mu _{n,m}\)  and \(b^{(2)} \sim \mu ^{\text{edit}}_{n, m} = \mu _{n+k_2,m}\) .
[Proof of claim:qsgdcoupling1]
\(b^{(1)} \sim \mu _{n,m}\)  follows trivially by construction. For the other part, for any set \(E\)  of \(m\)  indices arising from the coupling construction, let \(E_1\)  be the set of \(m-m_3\)  points from \(S \cup S_2^{\prime }\)  and \(E_2\)  be the set of \(m_3\)  points from \(S_3^{\prime }\) . Since these \(m_3\)  points of \(E_2\)  need to be selected when sampling \(b\) , the probability of sampling these points is \(\frac{{n \atopwithdelims ()m-m_3}}{{n+k_2 \atopwithdelims ()m}}\) , where the numerator denotes the number of ways to sample from \(S^{\prime } \backslash S_3^{\prime }\) . For the points in \(E_1\) , these come from \(b^{(1)}\)  and replacement using \(S_2^{\prime }\)  (which is a deterministic operation).
Hence, probability of \(E_1\)  is \(\frac{{n-(m-m_3) \atopwithdelims ()m_3}}{{n \atopwithdelims ()m}}\) , where the numerator denotes the number of ways to sample rest of elements not in \(E_1\)  when sampling \(b^{(1)}\) . Finally, we need to consider the re-sampling step i.e sampling \(b^{(2)}_1\)  - note that the draw of \(E_1\)  and \(E_2\)  fixes the set produced by this re-sampling, and thus its probability is \(\frac{1}{{m \atopwithdelims ()m_3}}\) . This gives us
\(\mu ^{\text{edit}}_{n, m}(E) &= \frac{{n \atopwithdelims ()m-m_3}}{{n+k_2 \atopwithdelims ()m}} \cdot \frac{{n-(m-m_3) \atopwithdelims ()m_3}}{{n \atopwithdelims ()m}} \cdot \frac{1}{{m \atopwithdelims ()m_3}} \\&= \frac{1}{{n+k_2 \atopwithdelims ()m}} \cdot \frac{n! (n-(m-m_3))! m! (n-m)! m_3! (m-m_3)!}{(m-m_3)! (n-(m-m_3))! m_3! (n-m)! n! m!} = \frac{1}{{n+k_2 \atopwithdelims ()m}} = \mu _{n+k_2,m}(E)\) 
In the second case, \(S\)  has more samples than \(S^{\prime }\)  - let number of samples in \(S\)  be \(n\) , and in \(S^{\prime }\)  be \(n-k_2\)  and there \(k_1\)  samples in \(S^{\prime }\)  not in \(S\) . As before we order the sets as: let \(S_1=S_1^{\prime }\)  be the \(n-k_2-k_1\)  samples which are the same in both \(S\)  and \(S^{\prime }\) . Let \(S_2\)  be the next \(k_1\)  samples in \(S\) , which correspond to \(S_2^{\prime }\) , the rest of \(k_1\)  samples in \(S^{\prime }\)  - the mapping from \(S_2\)  to \(S_2^{\prime }\)  being \(\iota \) . Finally let \(S_3\)  be the rest of \(k_2\)  samples in \(S\) .
We first sample \(b^{(1)} \sim \mu _{n,m}\)  (during training).
Let \(b^{(1)}_2\)  be the \(m_2\)  indices in \(S_2\) : replace these by the corresponding indices in \(S_2^{\prime }\)  i.e. \(\iota (b^{(1)} \cap S_2)\) .
Let \(b^{(1)}_3\)  denote the sub-sampled indices which are in the last \(k_2\)  indices of \(S\) , and let \(m_3 = {b^{(1)}_3}\) . We re-sample \(m_3\)  indices as \(b = \text{Unif}((S^{\prime } \backslash \iota (b^{(1)} \cap S_2) \backslash b^{(1)}),m_3)\) .
Finally, define \(b^{(2)} = {b^{(1)} \backslash {b^{(1)} \cap S_2} \cup {\iota (b^{(1)} \cap S_2)}}\backslash b^{(1)}_3 \cup b\) . Let the distribution of \(b^{(2)}\)  produced in the above way be denoted as \(\mu ^{\text{edit}}_{n, m}\) .
We now show that \((b^{(1)}, b^{(2)})\)  is indeed a coupling of \(\mu _{n,m}\)  and \(\mu _{n-k_2,m}\) .
With the construction described above, we have that \(b^{(1)} \sim \mu _{n,m}\)  and \(b^{(2)} \sim \mu ^{\text{edit}}_{n, m} = \mu _{n-k_2,m}\) .

[Proof of claim:qsgdcoupling2]
\(b^{(1)} \sim \mu _{n,m}\)  follows trivially by construction. For the other part, let \(E\)  be a set of \(m\)  indices from \([n-k_2]\) .
Note that any number of points in \(E\)  can arise due to re-sampling (i.e. when sampling \(b^{\prime }\) ), hence we need to consider all such possibilities - let \(m_3\)  be the number of indices in \(E\)  produced via re-sampling. Fixing one of \({m \atopwithdelims ()m_3}\)  combinations, the probability that it was re-sampled is \(\frac{1}{{n-k_2-(m-m_3) \atopwithdelims ()m_3}}\) . From the rule of sum, the probability that any \(m_3\)  sized set was produced via re-sampling is \(\frac{{m \atopwithdelims ()m_3}}{{n-k_2-(m-m_3) \atopwithdelims ()m_3}}\) .
For each such set, it could arise from any of \(m_3\)  points from \(k_2\) , which gives us \({k_2 \atopwithdelims ()m_3}\)  possibilities. The probability of choosing any such set, when sampling \(b^{(1)}\) , is \(\frac{{k_2 \atopwithdelims ()m_3}}{{n \atopwithdelims ()m}}\) . We now combine these and apply the rule of sum on different choices of \(m_3\) , from 0 to \(m\) . We get,
\(\mu ^{\text{edit}}_{n, m} (E) &= \sum _{m_3=0}^m \frac{{m \atopwithdelims ()m_3}}{{n-k_2-(m-m_3) \atopwithdelims ()m_3}} \cdot \frac{{k_2 \atopwithdelims ()m_3}}{{n \atopwithdelims ()m}} \\& = \frac{1}{{n \atopwithdelims ()m}}\sum _{m_3=0}^m\frac{m!(n-k_2-m)!m_3! }{m_3!(m-m_3)! (n-k_2-m+m_3)!} {k_2 \atopwithdelims ()m_3}\\& = \frac{1}{{n \atopwithdelims ()m}}\sum _{m_3=0}^m \frac{m!(n-k_2-m)!}{(n-k_2)!} \cdot \frac{(n-k_2)!}{(m-m_3)! (n-K_2-(m-m_3))!}{k_2 \atopwithdelims ()m_3}\\& = \frac{{n-k_2 \atopwithdelims ()m}}{{n \atopwithdelims ()m}} \sum _{m_3=0}^m {n-k_2 \atopwithdelims ()m-m_3} {k_2 \atopwithdelims ()m_3} = \frac{{n-k_2 \atopwithdelims ()m}}{{n \atopwithdelims ()m}} \cdot {n \atopwithdelims ()m} ={n-k_2 \atopwithdelims ()m} = \mu _{n-k_2,m}\) 

where the third last equality follows from Vandermonde's identity.
[!htpb]
Batch unlearning for quantized-m-SGD

[1]
Edit request produces dataset \(S^{\prime }\)  of \(n+k_2\)  points, with \(k_1\)  deletions and \(k_1+k_2\)  insertions; let \(S_1, S_2\)  and \(S_1^{\prime }, S_2^{\prime }\)  and \(S_3^{\prime }\)  be partitions of \(S\)  and \(S^{\prime }\)  respectively, as defined in “Batch coupling"
\(t=1,2\ldots , T\) 
Load(\(\theta _t,_t, b_t, _t\) )
\(b \sim \text{Unif}(S^{\prime },m)\) 
\(m_3 = {{x \in S_3 \cap b}}\) 
\(b^{(2)}_1 \sim \text{Unif}(b_t,m_3)\) 
\(_t^{\prime } = _t - \frac{1}{m}{\sum \limits _{j \in b_t \cap S_2}\nabla f(_t,_j) + \sum \limits _{j \in \iota (b_t \cap S_2)} \nabla f(_t,_j) - \sum \limits _{j \in b^{(2)}_1}\nabla f(_t,_j) + \sum \limits _{j \in S_3 \cap b}\nabla f(_t,_j)} \) 
\(b^{(2)}_t = {b_t \backslash {b_t \cap S_2} \cup {\iota (b_t \cap S_2)}}\backslash b^{(2)}_1 \cup {S_3^{\prime } \cap b}\) 
Save(\( _t^{\prime }, b^{(2)}_t\) )
\(Q_{\theta _t}{_t} \ne Q_{\theta _t}{_t^{\prime }}\) 
quantized-m-SGD // Recompute on current dataset
We now state the main result about unlearning.
(Algorithm REF , Algorithm REF ) satisfies exact batch unlearning.
Moreover, for \(k\)  batch edit requests, where the \(i^{\text{th}}\)  request comprises of \(k^i_1\)  deletions and \(k^i_1+k^i_2\)  insertions, or \(k^i_1\)  insertions and \(k^i_1+k^i_2\)  deletions, Algorithm REF  recomputes with probability at most \(2\sum _{i=1}^k(k^i_1+k^i_2)\rho \) .
[Proof of prop:unlearn-qsgd]
We consider one batch edit request of \(k_1\)  deletions and \(k_1+k_2\)  insertions (case 1) and \(k_1\)  insertions and \(k_1+k_2\)  deletions (case 2).
We have that applications of Claims REF  and REF  give us that mini-batches are transported, for cases 1 and 2 respectively.
Moreover, since we consider a trivial coupling of quantization randomness, we can consider it part of the (randomized) algorithmic map. Therefore, as in the proof of prop:unlearn-subsample-GD, transportation of mini-batches suffices to give us that alg:unlearnqsgd satisfies exact unlearning.
Repeated application of the above generalizes it to arbitrary \(k\)  edits.
We now proceed to bound the probability of recompute directly for a batch edit request.
For a fixed model \(\) , and a fixed iteration, we fix
the mini-batches \((b^{(1)},b^{(2)})\)  such that \(b^{(1)}\)  and \(b^{(2)}\)  differ by \(j\)  indices.
From claim:quantize, we have
\(\mathbb {P}_{\theta ,(b^{(1)},b^{(2)})}\left[Q_\theta (_{b^{(1)}}^S()) \ne Q_\theta (_{b^{(2)}}^{S^{\prime }}()) | (b^{(1)},b^{(2)}) \text{ such that }b^{(1)}, b^{(2)}\text{ differ in }j \text{ indices} \right]\le \frac{4Gdj}{m\tau }\) 
We now integrate over the conditioning event. To do this, we need to compute the probability of the event that sampling \((b^{(1)},b^{(2)})\)  generates \(j\)  differing indices - denote this as \(p(j)\) .
Since we have two case for coupling constructions, we consider each one by one. We first look at the second case:
from construction of the coupling, it is easy to verify that \(j\)  differing indices can be produced when, for any \(i\) , \(b^{(1)}\)  samples \(i\)  elements from the \(k_1\)  differing items and \(j-i\)  indices from the last \(k_2\)  indices, for any \(i\)  from 0 to \(j\) . Hence, by direct computation, we have
\(p(j) = \sum _{i=0}^j\frac{{k_1 \atopwithdelims ()i} {k_2 \atopwithdelims ()j-i} {n-(k_1+k_2) \atopwithdelims ()m-j}}{{n \atopwithdelims ()m}} = \frac{{n-(k_1+k_2) \atopwithdelims ()m-j}}{{n \atopwithdelims ()m}}\sum _{i=0}^j{k_1 \atopwithdelims ()i} {k_2 \atopwithdelims ()j-i} = \frac{{n-(k_1+k_2) \atopwithdelims ()m-j}{k_1+k_2 \atopwithdelims ()j}}{{n \atopwithdelims ()m}}\) 
where the last equality follows from Vandermonde's identity. Plugging this in the following, we have,
\(&\mathbb {P}_{\theta ,(b^{(1)},b^{(2)})}\left[Q_\theta (_{b^{(1)}}^S()) \ne Q_\theta (_{b^{(2)}}^{S^{\prime }}())\right] \\&= \sum _{j=0}^{k_1+k_2} p(j) \mathbb {P}_{\theta ,(b^{(1)},b^{(2)})}\left[Q_\theta (_{b^{(1)}}^S()) \ne Q_\theta (_{b_2}^{S^{\prime }}()) |b^{(1)}, b^{(2)}\text{ differ in }j \text{ indices} \right] \\&\le \sum _{j=0}^{k_1+k_2}\frac{{n-(k_1+k_2) \atopwithdelims ()m-j}{k_1+k_2 \atopwithdelims ()j}}{{n \atopwithdelims ()m}} \frac{4Gdj}{m\tau }= \frac{m(k_1+k_2)}{n}\frac{4Gd}{ m\tau } \\& = \frac{4Gd(k_1+k_2)}{n\tau } \le \frac{k\rho }{T}\) 

where the second equality is a consequence of Vandermonde's identity proved in lemma:tv-stability-qsgd (Base case \(k_2=0\) ) and the last inequality follows by plugging in \(\tau = \frac{4GdT}{\rho n}\) .
We now look at the first case (when \(S\)  is smaller than \(S^{\prime }\) ), which is slightly more involved. Let \(i_1\)  denote the number of indices in \(b^{(1)} \cap S_1\) ,
and let \(i_2\)  be the number of indices in \(b^{(2)} \cup S_3^{\prime }\) . Furthermore, since we resample \(i_2\)  indices from \(b^{(1)}\) , let \(i_3\)  be the number of indices from \(S_1\)  which are re-sampled. It can be verified that if \(b^{(1)}\)  and \(b^{(2)}\)  differ in \(j\)  indices, then we need to have \(i_1+i_3=j\) . This is because it can happen that both \(i^{(2)}\)  is large, but upon re-sampling, it chooses elements from \(k_1\) , which does not increase the number of different indices between \(b^{(1)}\)  and \(b^{(2)}\)  Also, note that by construction \(i_3\le i_2\le j\) . Hence the probability \(p(j)\) , by direct computation is,
\(p(j) &= \sum _{i_1,i_2,i_3=0, i_1+i_3=j, i_3\le i_2\le j}^j\frac{{k_1 \atopwithdelims ()i_1}{n-k_1\atopwithdelims ()m-i_1}}{{n \atopwithdelims ()m}}\cdot \frac{{k_2 \atopwithdelims ()i_2} {n \atopwithdelims ()m-i_2}}{{n+k_2 \atopwithdelims ()m}} \cdot \frac{{i_1 \atopwithdelims ()i_2-i_3}}{{m \atopwithdelims ()i_2}} \\& = \sum _{i_1=0}^j\sum _{i_2=j-i_1}^j\frac{{k_1 \atopwithdelims ()i_1}{n-k_1\atopwithdelims ()m-i_1}}{{n \atopwithdelims ()m}}\cdot \frac{{k_2 \atopwithdelims ()i_2} {n \atopwithdelims ()m-i_2}}{{n+k_2 \atopwithdelims ()m}} \cdot \frac{{i_1 \atopwithdelims ()i_2-(j-i_1)}}{{m \atopwithdelims ()i_2}}\) 

where in the second equality, we substituted \(i_3 = j-i_1\) . We now claim that \(p(j) = \frac{{n-k_1 \atopwithdelims ()m-j}{k_1+k_2 \atopwithdelims ()j}}{{n+k_2 \atopwithdelims ()m}}\) , which we will argue via a double counting argument. Note that it suffices to show that
\( \sum _{i_1=0}^j\sum _{i_2=j-i_1}^j{k_1 \atopwithdelims ()i_1}{n-k_1\atopwithdelims ()m-i_1}\cdot {k_2 \atopwithdelims ()i_2} {n \atopwithdelims ()m-i_2}\frac{{i_1 \atopwithdelims ()i_2-(j-i_1)}}{{m \atopwithdelims ()i_2}} = {n-k_1 \atopwithdelims ()m-j}{k_1+k_2 \atopwithdelims ()j}{n \atopwithdelims ()m}\) 
. Consider set \(A\)  of \(n+k_2\)  elements, composed of \(A_1\)  of \(n-k_1\) , \(A_2\)  of \(k_1\)  and \(A_3\)  of \(k_2\)  elements, and a \(B\)  of \(n\)  elements, composed of \(B_1\)  of \(n-k_1\)  and \(B_2\)  of \(k_1\)  elements. Note that the expression \({n-k_1 \atopwithdelims ()m-j}{k_1+k_2 \atopwithdelims ()j}{n \atopwithdelims ()m}\)  is the size of number of combinations of \(2m\)  elements, \(m\)  each from \(A\)  and \(B\)  such that the number of elements from \(A_2 \cup A_3\)  is \(j\) .
We will show that the other expression also counts this set, via basic combinatorial rules. For this, consider combinations of \(m\)  elements from \(B\)  \(A_3\)  and such that we have \(i_2\)  elements from \(A_3\)  and the rest \(m-i_2\)  from \(B\) . Also, consider combinations of \(m\)  elements from \(A_1 \cup A_2\)  which consists of \(i_1\)  elements from \(A_2\)  the rest from \(A_1\) . We now modify these as follows, out of \(m\)  elements from \(A\) , select \(i_2\)  elements and replace thse from elements from \(A_3\)  - not that if it turns out that out of \(i_2\)  selected, \(j-i_1\)  are from \(A_1\) , then the number of elements from \(A_2 \cup A_3\)  after replacement becomes exactly \(j\) . However, also note that for each such combination arising, there are \({m \atopwithdelims ()i_2}\)  combinations of samples from \(A_1\)  and \(A_2\) , which give the same final combination after replacement. Hence, we need to apply the rule of division, so as not to repeatedly count the same combination. Finally,
using the rule of sum to consider all possible values of \(i_1\)  and \(i_2\)  retrieves the expression \( \sum _{i_1=0}^j\sum _{i_2=j-i_1}^j{k_1 \atopwithdelims ()i_1}{n-k_1\atopwithdelims ()m-i_1}\cdot {k_2 \atopwithdelims ()i_2} {n \atopwithdelims ()m-i_2}\frac{{i_1 \atopwithdelims ()i_2-(j-i_1)}}{{m \atopwithdelims ()i_2}} = {n-k_1 \atopwithdelims ()m-j}{k_1+k_2 \atopwithdelims ()j}{n \atopwithdelims ()m}\)  and completes the argument.
We again plug in the above in the following expression to get,
\(&\mathbb {P}_{\theta ,(b^{(1)},b^{(2)})}\left[Q_\theta (_{b^{(1)}}^S()) \ne Q_\theta (_{b^{(2)}}^{S^{\prime }}())\right] \\&= \sum _{j=0}^{k_1+k_2} p(j) \mathbb {P}_{\theta ,(b^{(1)},b^{(2)})}\left[Q_\theta (_{b^{(1)}}^S()) \ne Q_\theta (_{b_2}^{S^{\prime }}()) |b^{(1)}, b^{(2)}\text{ differ in }j \text{ indices} \right] \\&\le \sum _{j=0}^{k_1+k_2}\frac{{n-k_1 \atopwithdelims ()m-j}{k_1+k_2 \atopwithdelims ()j}}{{n+k_2 \atopwithdelims ()m}} \frac{4Gdj}{m\tau }= \frac{m(k_1+k_2)}{n}\frac{4Gd}{ m\tau } \\& = \frac{4Gd(k_1+k_2)}{n\tau }\) 

where the second equality is again a consequence of Vandermonde's identity as in lemma:tv-stability-qsgd, and the last inequality follows by plugging in \(\tau = \frac{4GdT}{\rho n}\) . Finally, we condition on the iterates till iteration \(t\) , which gives us the conditional probability of the iterates differing at iteration \(t\)  is at most \(\frac{(k_1+k_2)\rho }{T}\) . Taking a union bound over all \(T\)  iterations
gives us that probability is at most \((k_1+k_2)\rho \) . Finally, we extend it to \(k\)  edit request, by using the fact, by assumption than the number of data points at any point in the stream is between \(\frac{n}{2}\)  and \(2n\) . This, with the result for one edit request, directly give us the probability to recompute is at most \(2\sum _{i=1}^k(k^i_1+k^i_2)\rho \) .
We now state and prove the main result.
Let \(f(\cdot , )\)  be an \(L\) -smooth \(G\) -Lipschitz convex function \(\forall \ \) .
For any \(\frac{1}{n} \le \rho < \infty \) , using alg:qsgd as the learning algorithm and alg:unlearnqsgd as its unlearning algorithm, then given a stream of batch edit requests,

Satisfies exact batch unlearning at every point in the stream.

At time \(i\)  in the stream of edit requests, outputs \(\hat{}_{S^i}\) , such that its excess empirical risk bounded as,
\({\hat{F}_S(\hat{}_{S^i}) - \hat{F}_{S}(_{S^i}^*)} \lesssim {\frac{\sqrt{L}GD^2 d^{3/2}}{\rho n}}^{2/3}\) 

For \(k\)  batch edit requests, where the \(i^{\text{th}}\)  request comprises of \(k^i_1\)  deletions and \(k^i_1+k^i_2\)  insertions, or \(k^i_1\)  insertions and \(k^i_1+k^i_2\)  deletions, the expected total unlearning runtime is
\(O(\max {\min {\rho ,1} \sum _{i=1}^k (k_1^1+k_2^i) \cdot \text{Training time},k})\)

[Proof of thm:main-thm-qsgd]
The first and the second claims follow from prop:unlearn-qsgd and prop:upper-bound-qsgd respectively combined with the assumption that the number of samples at every point in the stream is between \(\frac{n}{2}\)  and \(2n\) .
Finally, as in the the proof claim:runtimenoisyASGD for runtime noisy-m-A-SGD, we can use the same data-structures together with the fact the quantization operation takes \(O(d)\)  time, to get that the claimed runtime. These together finish the proof of thm:main-thm-qsgd.

Lower bounds on excess empirical risk
Give a convex function \(f(\cdot , )\) , we consider empirical risk minimization on a dataset of \(n\)  points.
We assume \(f(\cdot , )\)  is 1-Lipschitz for all \(\) , and diam\(() \le 1\) . This is only for simplification as the bounds scale naturally with these constants, as discussed in [15]}. We look at algorithms, which given two datasets \(S\)  and \(S^{\prime }\)  of size \(n\)  differing by one point, disagree only on a set of measure at most an \(\rho \) .
We have from the optimal transport connection that this requirement is equivalent to the total variation distance being at most \(\rho \) . We want to understand then what is the lower bound on excess empirical risk:
\(\underset{\Delta (S,S^{\prime })=1}{\sup } \mathbb {P}\left[(S^{\prime }) \ne (S)\right] \le \rho &\iff \underset{\Delta (S,S^{\prime })=1}{\sup }\text{TV}((S),(S^{\prime })) \le \rho \\&\Rightarrow {\left[\text{excess empirical risk}\right]} \ge \alpha (\rho ,n,d)\) 
We focus on proving the implication. [15]} gave lower bounds on accuracy for DP algorithms by providing a reduction to computing mean of the dataset. We present and give the proof of the reduction, adapted to our context, for completeness. The reduction is that if we have a TV-stable algorithm for empirical risk minimization for a particular \(f\)  with some accuracy, then we have a TV-stable algorithm for mean computation problem with certain accuracy.
We will look at mean computation problem over datasets with norm of the mean being \(\Theta (M)\) , for some given \(M\) . Let \(\mu (S) = \frac{1}{n}\sum _{j=1}^n_i\)  denote the mean of dataset \(S={_1,_2,\cdots _n}\) .
Let the optimal accuracy of such a mean computation problem be denoted as follows:
\(\alpha _{\text{mean}}^2(n,\rho , d,M):=\min _{: \rho \text{-TV-stable}} \max _{\begin{array}{c}S=\lbrace _i\rbrace _{i\in [n]}: {_i}\le 1,\\ M/2\le {\mu (S)}\le 2M\end{array}} _{(S)-\frac{1}{n}\sum _{i=1}^n_i}^2\) 
For any \(\rho \) -TV stable algorithm \(\) , there exists a 1-Lipschitz convex function \(f\) , a constraint set \(\)  with diameter(\()\le 1\)  and a dataset \(S\)  of \(n\)  data point such that
\(\hat{F}_S((S)) - \hat{F}_S(^*) \ge \max _M{\frac{\alpha _{\text{mean}}^2(n,\rho , d,M)}{2M}}\) 
[Proof of Proposition ]
We follow the proof in [15]}.
Consider dataset \(S={_1,_2,\ldots , _n}\) , \(_i \in {-\frac{1}{\sqrt{d}},\frac{1}{\sqrt{d}}}^d\)  - the dataset is therefore constrained to lie in the unit Euclidean ball. Consider the following function \(f(,) = -{}{}\)  with the constraint set \(\)  being the unit Euclidean ball. It is easy to see that \(f(\cdot ,)\)  is 1-Lipschitz for all \(\) . The empirical risk becomes \(\hat{F}_S() = -{}{\frac{1}{n}\sum _{i=1}^n_i}\) , the minimum of which over the unit ball is \(_S^* = \frac{\frac{1}{n}\sum _{i=1}^n_i}{{\frac{1}{n}\sum _{i=1}^n_i}}\) .
Given an algorithm \(\)  for empirical risk minimization, let the reduced mean estimate be \(\hat{\mu }(S) = {\frac{1}{n}\sum _{i=1}^n_i} (S)\) . The accuracy (mean-squared error) of \(\hat{\mu }\)  is,
\({\hat{\mu }(S) - \mu (S)}^2 &= {{\frac{1}{n}\sum _{i=1}^n_i} (S) - \frac{1}{n}\sum _{i=1}^n _i}^2= {{\frac{1}{n}\sum _{i=1}^n_i} {(S) - \frac{\frac{1}{n}\sum _{i=1}^n _i}{{\frac{1}{n}\sum _{i=1}^n_i}}}}^2 \\& = {\frac{1}{n}\sum _{i=1}^n_i}^2 {(S) - _S^*}^2\le {\frac{1}{n}\sum _{i=1}^n_i} 2{\hat{F}_S((S)) - \hat{F}_S(^*)}\) 

where the last inequality follows using the following computation, wherein we use the fact all data point are in the unit ball.
\({\frac{1}{n}\sum _{i=1}^n_i}^2 {(S) - _S^*}^2 &\le {\frac{1}{n}\sum _{i=1}^n_i} 2{1 - {(S)}{_S^*}} \\&= 2{{\frac{1}{n}\sum _{i=1}^n_i} - {\frac{1}{n}\sum _{i=1}^n_i}{(S)}{\frac{\frac{1}{n}\sum _{i=1}^n_i}{{\frac{1}{n}\sum _{i=1}^n_i}}}} \\& = 2{{^*_S}{\frac{1}{n}\sum _{i=1}^n_i} - {(S)}{\frac{1}{n}\sum _{i=1}^n_i}} \\&= 2(\hat{F}_S((S)) - \hat{F}_S(^*)\) 

We therefore get,
\(\hat{F}_S((S)) - \hat{F}_S(^*) \ge \frac{1}{2{\frac{1}{n}\sum _{i=1}^n_i}}{\hat{\mu }(S) - \mu (S)}^2\) 
There are two things left to show: a bound on \(\frac{1}{2{\frac{1}{n}\sum _{i=1}^n_i}}{\hat{\mu }(S) - \mu (S)}^2\)  and show that
the reduced algorithm \(\hat{\mu }(S)\)  is also \(\rho \) -TV stable.
We proceed with the latter:
note that \(\hat{\mu }(S) ={\frac{1}{n}\sum _{i=1}^n_i} (S)\) . However the term \({\frac{1}{n}\sum _{i=1}^n_i}\)  depends on the dataset, and even if, for a neighbouring dataset \(S^{\prime }\) , \((S)\)  and \((S^{\prime })\)  are \(\rho \) -close in total variation, this data dependent scaling can potentially increase the distance. However, if instead we define \(\hat{\mu }(S) =M (S)\) , where \(M\)  is a constant, then
it is indeed \(\rho \)  TV stable. Moreover, for reasonable values
of \(M\) , the there exists dataset for which \({\frac{1}{n}\sum _{i=1}^n_i} = \Theta (M)\) .
Finally, note that by definition, \({\hat{\mu }(S) - \mu (S)}^2 \ge \alpha _{\text{mean}}^2(n,\rho , d,M)\) . Taking a max over all \(M\)  gives us the desired statement.
Lower bound for mean computation
In this section, we look at the problem of mean computation with TV stability constraint.
Note that to establish lower bounds on excess empirical risk, we need to look at mean computation over data sets with means between \(M/2\)  and \(2M\) , for a given \(M\) . However, we will see the mean computation even over the unit ball has same accuracy convex ERM. We will therefore establish lower bounds for the general mean computation problem, but the construction will use datasets with means \(\Theta (M)\)  for certain values of \(M\) .
Given a dataset \(S={_1,_2,\ldots , _n}\) , where \(_i \in {-\frac{1}{\sqrt{d}},\frac{1}{\sqrt{d}}}^d\)  for all \(i,\)  the task is to compute the mean \(\mu (S) = \frac{1}{n}\sum _{j=1}_j\) , while ensuring that the procedure is \(\rho \) -TV-stable.
This task is often considered in the differential privacy literature, however with the data points being \(_j \in {0,1}^d\) . The mean computation task then corresponds to releasing all one-way marginals of the database. Since we want to consider data points which lie inside the Euclidean ball, we therefore scale it accordingly. Given an algorithm \((S)\) , the accuracy is defined as mean-squared error: \(\alpha ^2 := \alpha ^2_{\text{mean}}(n,\rho ,d) ={{(S)-\mu (S)}^2}\)  where the expectation is over the randomization of the algorithm.
We first describe two algorithms for this problem and give upper bounds.
Subsample-mean:
Consider an algorithm which sub-samples a \(\rho \) -fraction of the dataset and outputs the mean on it.


The Subsample-mean procedure satisfies \(\rho \) -TV-stability and has accuracy \(\alpha ^2 \le O{\frac{1}{\rho n}}\) .

[Proof of claim:meancomputationsubsampleupperbound]
The \(\rho \) -TV stability claim follows since TV distance is witnessed by the event that a differing sample is sub-sampled, which happens with probability \(\frac{\rho n}{n} =\rho \) .
The proof of accuracy follows from the proof of prop:upperboundsubsample, wherein we computed the gradient on uniformly sub-sampled \(m\)  out of \(n\)  points - we showed that the mean on sub-sampled points is an unbiased estimate of the average gradient. Furthermore, since the gradients were bounded as well, the expected accuracy of mean computation is the same as the variance of gradient computation, which we derived to be \(O{\frac{1}{m}} = O{\frac{1}{\rho n}}\) .

Noisy-mean:
The algorithm computes the mean and adds \(N(0,\sigma ^2 I)\)  noise to it with \(\sigma ^2 = \frac{C}{n^2\rho ^2}\) , where \(C\)  is an appropriate universal constant.
The Noisy-mean procedure satisfies \(\rho \) -TV-stability and has accuracy \(\alpha ^2 \le O{\frac{d}{(\rho n)^2}}\) 

[Proof of claim:meancomputationnoisyupperbound]
Since the difference in means of two datasets, in norm, is at most \(\frac{2}{n}\) , the outputs are two multivariate Gaussians with variance \(\sigma ^2\)  and means separated by \(\frac{2}{n}\) . From [1]}, the total variation distance between such Gaussian sis at most \(O{\frac{1}{n\sigma }} = \rho \) . For the accuracy, we have \(\alpha ^2 = {{\mu (S) + \xi - \mu (S)}^2} = {{\xi }^2} = d \sigma ^2 = O{\frac{d}{n^2\rho ^2}}\) .
If the above procedures are optimal, then we expect a lower bound of \(\alpha ^2 \gtrsim \min {\frac{1}{\rho n}, \frac{d}{(\rho n)^2}}\) . Equivalently, for a fixed accuracy \(\alpha \) , we expect a sample complexity lower bound of \(n\gtrsim \min {\frac{1}{\rho \alpha ^2},\frac{\sqrt{d}}{\rho \alpha }}\) .

Lower bound I
In this section, we give a \(\Omega {\frac{1}{\alpha \rho }}\)  lower bound on sample complexity.
The key ingredient is the following result, where the proof is based on a simple reduction argument.
Suppose there exists a \(\rho \) -TV-stable algorithm such that for any dataset of \(n\)  points, it achieves an accuracy of \( \alpha \) . Then there exists a \(0.1\) -TV stable algorithm such the for any dataset of size \(\lceil 100n \alpha \rho \rceil \) , it achieves a \(0.1\) -accuracy.
[Proof of prop:lowerbounddpreduction]
Let \(n^{\prime } =\lceil 100n \alpha \rho \rceil \) . Consider a dataset \(S^{\prime }\)  of \(n^{\prime }\)  points. We construct a dataset \(S\)  of \(n\)  points by concatenating \(K = \lceil 0.1/\rho \rceil \)  copies of \(S^{\prime }\)  followed by \(\lceil \frac{n-Kn^{\prime }}{2}\rceil \)  copies of a constant sample, all ones \((\frac{1}{\sqrt{d}},...,\frac{1}{\sqrt{d}})\)  and \(\lceil \frac{n-Kn^{\prime }}{2}\rceil \)  copies of a constant sample, all ones \((-\frac{1}{\sqrt{d}},...,-\frac{1}{\sqrt{d}})\) .
Consider the algorithm wherein we compute the stable-mean on \(D^{\prime }\)  by \(^{\prime }\) , defined as computing the stable-mean on \(D\)  using \(\)  and adjusting:
\(^{\prime }(D^{\prime }) = \frac{n}{Kn^{\prime }}(D)\) 
Let \(\tilde{S}^{\prime }\)  be a neighbouring dataset of \(S^{\prime }\) .
By construction, note that \(S\)  and \(\tilde{S}\)  differ by \(K\)  samples. Furthermore, since the algorithm \(\)  on \(D\)  is \(\rho \) -TV stable, on \(K\) -neighbouring datasets, it is \(K\rho \)  = \(0.1\) -TV. This establish the stability part of claim.
The accuracy, by direct computation is \({{^{\prime }(D^{\prime })-\mu (D^{\prime })}^2} = \frac{n^2}{K^2n^{\prime 2}}{{(D)-\mu (D)}^2} \le (0.1)^2\) .
For the \(d\) -dimensional mean computation problem over the Euclidean ball, there exists a dataset \(S\)  of \(n\)  samples with mean \({\mu (S)}=\Theta {\frac{1}{\rho n}}\)  such that
the accuracy of any \(\rho \)  TV stable algorithm is \(\alpha \ge \Omega {\frac{1}{\rho n}}\) .

[Proof of thm:lowerbound11]
Even for accuracy \(\alpha _0= 0.1\)  accuracy and \(\rho _0=0.1\)  stability, we need at least one sample.
Hence, using prop:lowerboundreduction, we get that sample complexity is \(n\ge \Omega {\frac{1}{\rho \alpha }}\) , which equivalently gives the claimed accuracy lower bound.
Note that for this one-sample dataset \(S^{\prime }\) , \({\mu (S^{\prime })} =1\) .
Finally, from the reduction in prop:lowerboundreduction, the mean of dataset \(S\)  becomes \({\mu (S)} = \lceil \frac{0.1}{\rho n}\rceil \) , which finishes the proof.

Lower bound II
In this section, we will prove the \(\Omega {\frac{1}{\alpha ^2 \rho }}\)  lower bound. We first introduce a technical assumption.
For any dataset \(S\) , we assume that the probability distribution \((S)\)  is defined over the unit Euclidean ball, is absolutely continuous with respect to the uniform measure (in the unit Euclidean ball) and its
probability density function, with respect to the uniform measure, is bounded by \(K\)  in absolute value.
As a remark, the above assumption can also be stated with respect to the Lebesgue measure, but then we would get a scaling of \(\frac{\pi ^{d/2}}{\Gamma {1+\frac{d}{2}}}\) , which is the Lebesgue volume of the \(B_d(0,1)\) , to some of our terms. In order to simplify, we therefore use the uniform measure.
Let \(n\ge 72, \alpha \le \frac{1}{4}\)  and \(\frac{1}{n} \le \rho \le \frac{1}{4}\) .
Let \(\)  be any \(\rho \) -TV-stable algorithm satisfying Assumption REF  with \(K\le 2^d\) .
For large enough dimension \(d\) , there exists a dataset \(S\)  of \(n\)  points with \({\mu (S)} = \Theta {\frac{1}{\sqrt{\rho n}}}\)  such that accuracy is lower bounded as \(\alpha \ge \Omega {\frac{1}{\sqrt{\rho n}}}\) .
[Proof of thm:lowerbound2]
We will prove the result by contradiction.
Let "\(\text{Vol}\) " of a set refer to its volume with respect to the uniform measure on the unit ball. Consider the following high-dimensional setup.
Consider a dataset \(S\)  (or \(S^0\) ) which mean \(\mu (S)\)  such that \({\mu (S)}=\Theta {\frac{1}{\sqrt{\rho n}}}\) . It is easy to construct such datasets by considering points such that sum of \(n-\lceil \sqrt{\frac{n}{\rho }}\rceil \)  points is 0 and the rest of points is the same point repeated - this uses the assumption that \(\rho \ge \frac{1}{n}\) .
Now consider neighbouring datasets \(S^i\) 's, \(i\in [n]\)  such that the means of \(S^i\) 's are all \(\frac{1}{n}\)  far from that of \(S\) , in norm.
We also need that the means of any two datasets \({\mu (S^i)-\mu (S^j)} \ge \frac{1}{2n}\)  for \(i,j=0\)  to \(n\)  and \(i\ne j\) . It is easy to see the existence of such datasets, by considering the means of \(S^i\) 's in near orthogonal directions to that of \(S\) , which is possible when \(d\)  is large enough.
Suppose the algorithm \(\)  has expected error \(\alpha ^2\)  i.e. \({{(S)-\mu (S)}}^2 \le \alpha ^2\)  with \( 72 \le n \le \frac{1}{\alpha ^2\rho }\) .
Consider \(B_d{\mu (S),\frac{1}{K^{1/d}}}\) , the \(d\)  dimensional Euclidean ball centered at \(\mu (S)\)  of radius \(\frac{1}{K^{1/d}}\) .
From Markov's inequality, we have that \(¶{(S)\notin B_d{\mu (S),\frac{1}{K^{1/d}}}} \le ¶{(S)\notin B_d{\mu (S),\frac{1}{2}}} = ¶{{(S)-\mu (S)}^2 \ge \frac{1}{4}}\le \frac{{{(S)-\mu (S)}}^2}{(1/4)} \le 4\alpha ^2\) , where in the first inequality, we used the assumption \(K\le 2^d\) . Therefore, we have \(¶{(S) \in B_d{\mu (S),\frac{1}{K^{1/d}}}} \ge 1-4\alpha ^2\) .
We now setup some additional notation. Let \(A_i\)  denote the set \(B_d{\mu (S^i),\frac{1}{K^{1/d}}} \backslash { \cup _{j=0,j \ne i}^n B_d{\mu (S^j),\frac{1}{K^{1/d}}}}\)  i.e. the region in the ball \(B_d(\mu (S^i)\)  which is not contained in any of the other balls.
Let \(B_{ij}\)  denote the region of intersection between \(B_d{\mu (S^i),\frac{1}{K^{1/d}}}\)  and \(B_d{\mu (S^j),\frac{1}{K^{1/d}}}\)  where \(i\ne j\)  and \(i\)  and \(j\)  go from 0 to \(n\) . Note that set \(B_{ij}\)  is constituted of two spherical caps. By construction the centers of the intersecting spheres are at least \(\frac{1}{2n}\)  apart.
To study the properties of such a set, we define cap as the region in a \(d\)  dimensional sphere of radius \(\frac{1}{K^{1/d}}\)  which intersects with another sphere of the same radius but with centers being apart by \(1/2n\) .
From known results [19]}, the volume of cap is asymptotic to \(\text{Vol}{B_d{0,\frac{1}{K^{1/d}}}}{1-\Phi {\frac{\sqrt{d}}{2n}}}\)  as \(d \rightarrow \infty \) 
where \(\Phi \)  is the cumulative distribution function of a standard normal random variable. We therefore have that \( \lim _{d \rightarrow \infty } \text{Vol(\textbf {cap})} \sim \lim _{d \rightarrow \infty }\text{Vol}{B_d{0,\frac{1}{K^{1/d}}}}{1-\Phi {\frac{\sqrt{d}}{2n}}} = 0\) .
Furthermore, using Assumption REF , we have \(¶{(S) \in \text{cap}} \le K \text{Vol(\textbf {cap})}\sim K\text{Vol}{B_d{0,\frac{1}{K^{1/d}}}}{1-\Phi {\frac{\sqrt{d}}{n}}} \le K {\frac{1}{K^{1/d}}}^d{1-\Phi {\frac{\sqrt{d}}{n}}} = {1-\Phi {\frac{\sqrt{d}}{n}}}\)  as \(d\rightarrow \infty \) .
Since \(\Phi (t) = \mathbb {P}_{g \sim (0,1)}[g \le t]\) , we have that \(1 - \Phi (t) = \mathbb {P}_{g \sim (0,1)}[g > t] \le \frac{e^{-t^2/2}}{\sqrt{2\pi }t} \)  where the last inequality follows from standard bounds on tails of normal distribution (See Proposition 2.1.2 in [20]}). Therefore, we have \(¶{(S) \in \textbf {cap}} \lesssim \frac{ne^{-d/4n^2}}{\sqrt{d}}\) . For constant \(\epsilon >0\) , choosing \(d \gtrsim 4n^2\ln {{n/\epsilon }}\)  ensures that \(¶{(S) \in \textbf {cap}} \le \frac{\epsilon }{2n}\)  for large enough \(n\)  (to be specified later). Since \(B_{ij}\)  is made up of two conjoined caps, this gives us that for any \(j=0\)  to \(n\)  and \(i\ne j\) , we have that \(¶{(S^i)\in B_{ij}} \le \frac{\epsilon }{n}\) .
Finally, we look at \(A_i\) 's by removing the mass of all \(B_{ij}\) 's, and using a union bound, we get that \(¶{(S^i) \in A_i} = ¶{(S^i) \in B_d{\mu (S^i),\frac{1}{K^{1/d}}}} - ¶{ (S^i) \in \cup _{j=0,j\ne i}^n B_{ij}} \ge 1-4\alpha ^2 - \sum _{j=0,j\ne i}^{n} ¶{(S^i) \in B_{ij}} \ge 1-4\alpha ^2 -\epsilon \ge \frac{1}{2}-4\alpha ^2 \)  where the last inequality holds for \(\epsilon \le \frac{1}{2}\) . We now evaluate how large \(n\)  we need for this regime of \(\epsilon \) : recall that we set \(d\gtrsim 4n^2\ln {{n/\epsilon }}\) , this gives
\(\frac{ne^{-d/4n^2}}{\sqrt{d}} \le \frac{\epsilon }{2n\sqrt{ \ln {{n/\epsilon }}}}\) . We want the right hand side to be at most \(\frac{\epsilon }{2n}\)  for \(\epsilon \le 1/2\) . Plugging in this worst-case value of \(\epsilon \) , we get the condition \(\ln {(2n)} \ge 1\)  which holds for any \(n\ge 1.4\)  and therefore is valid by our assumption of \(n\) .
We now use the fact that \(A_i^{\prime }s\)  are disjoint by construction. Therefore the total measure of \((S)\)  on union of \(A_i^{\prime }s\)  is at most 1 i.e \( ¶{(S) \in \cup _{i=1}^n A_i} = \sum _{i=1}^n ¶{(S) \in A_i} \le 1\) . Furthermore, since \((S)\)  is \(\rho \) -TV stable, we have that \(¶{(S) \in A_i} \ge ¶{(S^i) \in A_i} - \rho \) .
Combining this with the previous analysis which gives a lower bound on \(¶{(S^i) \in A_i}\)  yields
\(n{\frac{1}{2}-4\alpha ^2-\rho } \le \sum _{i=1}^n¶{(S^i) \in A_i} - \rho \le \sum _{i=1}^n ¶{(S) \in A_i} \le 1\) 
We now proceed in two cases:
Case 1: Suppose \(72 \le n \le \frac{17}{4\alpha ^2}\) . The latter condition gives us that \(4\alpha ^2 \le \frac{17}{n}\) . Using eqn:lb2 gives us \(n(1/2-4\alpha ^2-\rho )\le 1 \iff 4\alpha ^2 \ge \frac{1}{2}-\frac{1}{n}-\rho \) . Upper bounding \(4\alpha ^2\)  by \(\frac{17}{n}\)  gives us that \(\frac{18}{n} \ge \frac{1}{2} - \rho \iff n\le \frac{18}{(1/2-\rho )} \le 72\)  where in the last inequality we used \(\rho \le 1/4\) . This gives us a contradiction.
Case 2: Suppose \( \frac{17}{4\alpha ^2} \le n \) . We again start with eqn:lb2 which gives us \(n(1/2-4\alpha ^2-\rho )\le 1 \iff \rho \ge \frac{1}{2}-\frac{1}{n}-4\alpha ^2\) . We want to prove the right hand side is at least \(\frac{1}{n \alpha ^2}\) , which would give us that \(n\ge \frac{1}{\rho \alpha ^2}\) . Suppose this is not true i.e. \(\frac{1}{2} - \frac{1}{n} -4\alpha ^2 \le \frac{1}{n \alpha ^2} \iff \frac{n-2-8\alpha ^2n}{2n} \le \frac{1}{n\alpha ^2} \iff \alpha ^2n(1-8\alpha ^2) \le 2(1+\alpha ^2) \iff n \le \frac{2(1+\alpha ^2)}{\alpha ^2(1-8\alpha ^2)}\) . Finally using the fact that \(\alpha \le \frac{1}{4}\)  gives that \(n \le \frac{2(1+1/16)}{\alpha ^2(1-8/16)} \le \frac{17}{4\alpha ^2}\)  which yields a contradiction.
Hence, we see that with \(n>72\)  samples and accuracy \(\alpha ^2\) , we have established that \(n\ge \frac{1}{\rho \alpha ^2}\)  and so \(\alpha \ge \frac{1}{\sqrt{\rho n}}\) .

Excess population risk bounds
The goal in machine learning is (population) risk minimization.
The population risk of \(\) , denoted by \(F()\)  is defined as \(F() := {\sim }{f(,)}\) , where \(\)  is an unknown probability distribution over data points. Analogously, given an output of algorithm \(\)  on dataset \(S ={_i}_i\)  where \(_i \sim \)  i.i.d., denoted as \((S)\) , we will give guarantees on the expected excess population risk, defined as \({F((S))- F(^*)}\) , where \(^*\)  is the population risk minimizer: \(^* \in \arg \min _{\in } F()\) , and the expectation is taken with respect to randomness in algorithm \(\)  as well as sampling \(S\) .
Upper bounds
In this section, we will bound the expected excess population risk appealing to connections between algorithmic stability and generalization [21]}. We first define uniform stability.
[Uniform stability]
Let \(:^n \rightarrow \)  be an algorithm and \((S)\)  denotes its output on dataset \(S\) . We say that \(\)  is \(\epsilon _{\text{stable}}(n)\) -uniformly stable if for any datasets \(S\)  and \(S^{\prime }\)  of \(n\)  points such that they differ by one data point (i.e. \(\Delta (S,S^{\prime })= 2\) ), we have \(\sup _{z\in }\mathbb {E}_{}\left[f((S),z)- f((S^{\prime }),z)\right]\le \epsilon _{\text{stable}}(n)\)
A classical result [21]} shows that expected excess population risk is at most uniform stability + expected excess empirical risk: i.e. any \(\in \) , we have
\(\mathbb {E}\left[F((S)) - F()\right] \le \epsilon _{\text{stable}}(n)+ \mathbb {E}\left[\hat{F}_S((S)) - \hat{F}_S() \right]\) 
[Upper bound]

There exists a \(\rho \)  TV stable algorithm, such that for any function \(f(\cdot , )\)  which is \(L\) -smooth \(G\) -Lipschitz convex \(\forall \ \)  and any dataset \(S\)  of \(n\)  points, it outputs \(\hat{}_S\)  which satisfies the following.
\({ F(\hat{}_S) - F(^*)} \lesssim \frac{GD}{\sqrt{n}}+ GD\min {\frac{1}{\sqrt{\rho n}}, \frac{\sqrt{d}}{\rho n}}\) 
[Proof of thm:pop-risk-upper-bound]
We use sub-sample-GD (alg:sub-sample-GD) and noisy-m-SGD (alg:noisy-m-sgd). From Lemma 3.2 in [23]}, we have that \(\epsilon _{\text{stable}}(n) \le \frac{G^2\eta T}{n}\) . From prop:upperboundsubsample, we set \(\eta = \min {\frac{1}{2L},\frac{D}{G\sqrt{T}}}\) , and \(T= \frac{DL\sqrt{\rho n}}{G}\) . We therefore have \(\epsilon _{\text{stable}}(n) \le \frac{G^2T}{2Ln} = \frac{GD\sqrt{\rho }}{\sqrt{\rho n}}\) . Using the excess empirical risk bound from prop:upperboundsubsample, and the fact that \(\rho \le 1\) , the excess population risk is bounded as,
\({ F(\hat{}_S) - F(^*)} \le \frac{GD\sqrt{\rho }}{\sqrt{n}} + \frac{GD}{\sqrt{\rho n}} \le \frac{GD}{\sqrt{n}}+\frac{GD}{\sqrt{\rho n}}\) 
For noisy-m-SGD, we need to balance the trade-offs more directly. In prop:upperboundnoisy-m-sgd, we arrived at that when using \(\eta \le \frac{1}{2L}\) , the expected excess empirical risk is bounded by \(\eta ^2 + \frac{D^2}{\eta T}\) . Using the uniform stability bound of \(\frac{G^2T \eta }{n}\) , the expected excess population risk is bounded as,
\({ F(\hat{}_S) - F(^*)} \le \frac{G^2T \eta }{n} + \eta ^2 + \frac{D^2}{\eta T} = \eta {\frac{G^2 T}{n} + ^2} + \frac{D^2}{\eta T}\) 

Define \(\tilde{G}^2= {\frac{G^2 T}{n} + ^2}\) , where \(^2 \lesssim G^2 + \sigma ^2d \lesssim G^2 + \frac{G^2 Td}{n^2 \rho }\) . Setting \(\eta = \min {\frac{1}{2L},\frac{D}{\tilde{G} \sqrt{T}}}\) , we get,
\({ F(\hat{}_S) - F(^*)} &\le \frac{LD^2}{T} + \frac{\tilde{G} D}{\sqrt{T}}\\&\lesssim \frac{LD^2}{T} +\frac{G\sqrt{T}D}{\sqrt{T}\sqrt{n}} + \frac{GD}{\sqrt{T}} + \frac{GD\sqrt{T}\sqrt{d}}{\rho n\sqrt{T}} \\&= \frac{LD^2}{T} + \frac{GD}{\sqrt{n}} + \frac{GD}{\sqrt{T}} +\frac{GD \sqrt{d}}{\rho n}\) 
Setting \(T = \max { \min {\sqrt{n},\frac{\rho n}{\sqrt{d}}}, \frac{LD}{G}\min {\sqrt{n},\frac{\rho n}{\sqrt{d}}}}\) , and combining the two results finishes the proof.

Lower bounds
In this section, we will prove a lower bound on excess population risk for any \(\rho \) -TV stable algorithm. As before, we will consider the Lipschitz constant \(G\)  and diameter \(D\)  to be both 1, as the bounds scale naturally with these constants. We first define the following quantity, which denotes the lower bound on expected excess empirical risk of \(\rho \) -TV-stable algorithm with \(n\)  points.
\(\hat{\alpha }(n,\rho ) := \inf _{: \rho \text{-TV-stable}}\sup _{S: {S}=n} \mathbb {E}_\hat{F}_S((S)) - \hat{F}_S(\hat{}_S)\) 
For the problem of stochastic convex optimization, there exists a data distribution \(\) , such that any \(\rho \) -TV-stable algorithm \(\)  incurs expected excess population risk, bounded as follows
\(\underset{S \sim ^n, }{\mathbb {E}}F((S)) - F(^*) \ge \max {\Omega {\frac{1}{\sqrt{n}}}, \hat{\alpha }(n,\rho ) }\) 
[Proof of Theorem REF ]
The \(\frac{1}{\sqrt{n}}\)  term follows directly since it is the lower bound for any algorithm, and so applies to \(\rho \) -TV stable algorithms as well.
We now focus on the second term \(\hat{\alpha }(n,\rho )\) 
The proof is based on a standard reduction argument: if there is \(\rho \) -TV stable algorithm, which with \(n\)  i.i.d. samples from any distribution, achieves an expected excess population risk less than \(\hat{\alpha }(n,\rho )\) , then there is an \(\rho \) -TV stable algorithm which achieves an expected excess empirical risk less than \(\hat{\alpha }(n,\rho )\)  on any dataset of \(n\)  samples. Since the latter contradicts the definition of \(\hat{\alpha }(n,\rho )\) , this gives us that the expected excess population risk is at least or equal to \(\hat{\alpha }(n,\rho )\) . We now focus on the proof of the reduction. Consider a dataset \(S\)  of \(n\)  points. Consider \(\)  as the following algorithm: sample \(n\)  i.i.d. samples from \(S\) , call this set \(\tilde{S}\) , and run some \(\rho \) -TV algorithm \(\tilde{}\)  on \(\tilde{S}\) . For a fixed \(\tilde{S}\) , from TV-stable property of \(\tilde{}\) , for any neighbouring dataset \(\tilde{S}^{\prime }\)  with one point differing, we have that \(\text{TV}(\tilde{}(\tilde{S}),\tilde{}(\tilde{S}^{\prime })) \le \rho \) .
Furthermore, using the group property of TV-stability, for any dataset \(\tilde{S}^{\prime }\) , we have \(\text{TV}(\tilde{}(\tilde{S}),\tilde{}(\tilde{S}^{\prime })) \le \Delta (S,S^{\prime })\rho \) . Using the maximal coupling characterization of total variation distance, we have that there exists a coupling \(\tilde{\pi }\)  of random variables \(\tilde{}(\tilde{S})\)  and \(\tilde{}(\tilde{S}^{\prime })\)  such that \(\text{TV}(\tilde{}(\tilde{S}),\tilde{}(\tilde{S}^{\prime })) = \mathbb {E}_{\tilde{\pi }}\mathbb {1}{\tilde{}(\tilde{S}) \ne \tilde{}(\tilde{S}^{\prime })}\) . We now show that the algorithm \(\)  is also \(\text{TV}\) -stable for dataset \(S\) .
Consider dataset \(S^{\prime }\)  of \(n\)  points which differs from \(S\)  in the first sample. We now generate \(\tilde{S}^{\prime }\)  by drawing \(n\)  i.i.d samples from \(S^{\prime }\) . For this, consider the following coupling: we draw \(n\)  i.i.d samples from \(S\) , call it \(\tilde{S}\) . For every draw of the first sample, replace it by the first sample of \(S^{\prime }\) , call it \(\tilde{S}^{\prime }\) . It is easy to check the \(\tilde{S}\)  and \(\tilde{S}^{\prime }\)  are i.i.d. samples from \(S\)  and \(S^{\prime }\)  respectively. We now proceed to show the \(\)  is \(\rho \) -TV stable. We will use the fact the total variation distance is at most the probability of disagreement under any coupling. The coupling \(\pi \)  we consider is that we first generate \(\tilde{S}\)  and \(\tilde{S}^{\prime }\)  using the aforementioned coupling, and then use the coupling \(\tilde{\pi }\)  which achieves total variation distance for worst-case fixed neighbouring datasets \(\tilde{S}\)  and \(\tilde{S}^{\prime }\) . We have,
\(\text{TV}((S),(S^{\prime })) &\le \mathbb {E}_{\pi }\mathbb {1} {(S) \ne (S^{\prime })} = \mathbb {E}_{\pi }\mathbb {1}{\tilde{}(\tilde{S}) \ne \tilde{}(\tilde{S}^{\prime })} \\& \le \mathbb {E}_{\pi } \sup _{\tilde{S}, \tilde{S}^{\prime }}\mathbb {1}{\tilde{}(\tilde{S}) \ne \tilde{}(\tilde{S}^{\prime })} \le \mathbb {E}_{\pi } \Delta (\tilde{S},\tilde{S}^{\prime })\rho = \rho \) 

where the last equality follows from direct computation of \(\Delta (\tilde{S},\tilde{S}^{\prime })\) : number of differing samples, under coupling \(\pi \) .
We now proceed to the accuracy guarantee. From straight-forward computation, the excess population risk, under the sampling of \(\tilde{S}\) , is \(\hat{F}_S(\tilde{}(\tilde{S})) - \hat{F}_S(^*_S)\)  - this is the excess empirical risk for dataset \(S\) . So if we have an upper bound on excess population risk using algorithm \(\tilde{}\) , we have an upper bound on excess empirical risk for dataset \(S\) . This completes the reduction argument and hence the proof.

Algorithms for approximate unlearning
We first define the notion of approximate unlearning based on differential privacy (DP).
[\((\epsilon ,\delta )\) -approximate-unlearning]
We say a procedure \((\mathbf {A},\mathbf {U})\)  satisfies \((\epsilon ,\delta )\) -approximate-unlearning unlearning  if for any
\(S, S^{\prime } \subset ^*\)  such that \(\Delta (S,S^{\prime })=1\)  and for any measurable event \(\in \text{Range}(\mathbf {U}) \cap \text{Range}(\mathbf {A})\) , with probability at least \(1-\delta \) ,
\(^{-\epsilon }¶{\mathbf {U}(\mathbf {A}(S),S^{\prime }\backslash S \cup S\backslash S^{\prime }) \in } \le ¶{\mathbf {A}(S^{\prime }) \in } \le e^{\epsilon } ¶{\mathbf {U}(\mathbf {A}(S),S^{\prime }\backslash S \cup S\backslash S^{\prime }) \in }\) 
We now define \((g,\epsilon ,\delta )\) -group differential privacy.
[\((g,\epsilon ,\delta )\) -group differential privacy]
An algorithm \(\mathcal {A}\)  satisfies \((\epsilon ,\delta )\) -differential privacy if for any two datasets \(S\)  and \(S^{\prime }\)  such that \(\Delta (S,S^{\prime })\le g\) , for any measurable event \(\in \text{Range}()\) , it satisfies
\(\mathbb {P}(\mathcal {A}(S)\in ) \le e^\epsilon \mathbb {P}(\mathcal {A}(S^{\prime })\in )+\delta \) 
[24]}
If an algorithm satisfies \((\epsilon ,\delta )\) -DP, then for any \(g \in \) , it satisfies \((g, g \epsilon , g e^{(g-1)\epsilon }\delta )\) -group differential privacy.
We now define privateCompute oracle which, basically is a differentially private solver for the said task.
[privateCompute\((S,\epsilon ,\delta \) ) oracle]
For a problem instance, given a dataset \(S\)  of \(n\)  points, and privacy parameters \(\epsilon \)  and \(\delta \) , a privateCompute oracle outputs a \((\epsilon ,\delta )\) -differentially private solution with accuracy \(\alpha _{\text{private}}(n,\epsilon ,\delta )\)
We now give a very simple algorithm (alg:interpolate) based on the observation above using privateCompute oracle calls.
[!htpb]
Approximate unlearning

[1]
\(\epsilon ,\delta ,\rho \) 
\(i \leftarrow 0\) 
\(\hat{}_S \leftarrow \)  PrivateCompute \({S,\rho \alpha , \beta _{\text{group}}{\left\lfloor \frac{1}{\rho }\right\rfloor ,\rho \alpha , \rho \beta }}\) 
// Observe \(k\)  edit requests
\(t=1,2,\ldots , k\) 
\(S_t \leftarrow \text{Update dataset}\) (edit request)
\(i+=1\) 
\(i=\left\lfloor \frac{1}{\rho } \right\rfloor \) 
\(\hat{}_{S_t} \leftarrow \)PrivateCompute\({S_t,\rho \alpha , \beta _{\text{group}}{\left\lfloor \frac{1}{\rho }\right\rfloor ,\rho \alpha , \rho \beta }}\) 
\(i \leftarrow 0\)
Given a set of \(n\)  data points to start with, and observing a stream of \(k\)  requests, at any time \(t\)  in the stream, the following hold about Algorithm :

It satisfies \((\epsilon ,\delta )\) -approximate unlearning.

The unlearning runtime for \(k\)  requests is at most \(2\rho k\)  privateCompute oracle calls.

The accuracy is at most \(\alpha _\text{private}{\frac{n}{2},\rho \epsilon ,\delta _{\text{group}}{\left\lfloor \frac{1}{\rho }\right\rfloor , \rho \epsilon , \rho \delta }}\) .

[Proof of claim:interpolate]
Consider a point \(t\)  in the stream, and let \(j\)  be such that \(j\left\lfloor \frac{1}{\rho }\right\rfloor \le t \le (j+1)\left\lfloor \frac{1}{\rho }\right\rfloor \) .
Since the algorithm uses privateCompute with parameters \(\rho \epsilon \)  and \(\delta _{\text{group}}{\left\lfloor \frac{1}{\rho }\right\rfloor , \rho \epsilon , \rho \delta }\) , it satisfies \((\rho \epsilon ,\delta _{\text{group}}{\left\lfloor \frac{1}{\rho }\right\rfloor , \rho \epsilon , \rho \delta }\)  differential privacy and hence \({\left\lfloor \frac{1}{\rho }\right\rfloor ,\alpha ,\delta }\) -group privacy.
Therefore, for any such \(t\) , since the number of requests after time \(j \left\lfloor \frac{1}{\rho }\right\rfloor \)  is less that or equal to \(\left\lfloor \frac{1}{\rho }\right\rfloor \) , this implies it satisfies \((\alpha ,\delta )\) -approximate unlearning.
For the second part of the claim, note that for \(k\)  updates, the number of times the algorithm calls privateCompute is \(\frac{k}{\left\lfloor \frac{1}{\rho }\right\rfloor }\) .
Note that \(\frac{1}{\rho } \ge 1\) , so if \(1 \le \frac{1}{\rho } < 2\) , then \(\left\lfloor \frac{1}{\rho } \right\rfloor =1\) , which gives that the update complexity is \(k \le 2k\) . However, if \(\frac{1}{\rho } \ge 2\) , we have that \(\frac{k}{\left\lfloor \frac{1}{\rho } \right\rfloor } \le \frac{k}{{\frac{1}{\rho }-1}} \le 2k\rho \) , which gives the update complexity is at most \(2 \rho k\)  in both cases.
For the third part of the claim, at time \(j\left\lfloor \frac{1}{\rho }\right\rfloor \le t \le (j+1)\left\lfloor \frac{1}{\rho }\right\rfloor \) ,
the private estimator is computed with \(n{j\left\lfloor \frac{1}{\rho }\right\rfloor )} \ge \frac{n}{2}\) , by assumption.
Moreover the privacy parameters of the algorithm are \(\rho \epsilon \)  and \(\delta _{\text{group}}{\left\lfloor \frac{1}{\rho }\right\rfloor , \rho \epsilon , \rho \delta }\)  which gives the claimed accuracy bound.
As an example, consider \(\rho = \frac{1}{\sqrt{k}}\) ,
we first do privatecompute with parameters \((\epsilon /\sqrt{k},,\delta _{\text{group}}(\epsilon /\sqrt{k},\delta /\sqrt{k},\lfloor \sqrt{k}\rfloor ))\) . Since after \(\lfloor \sqrt{k} \rfloor \)  edit requests, we would no longer satisfy the unlearning guarantee, so we now need to do privateCompute again. However note that we would only need to do privateCompute \(\sqrt{k}\)  times which gives the update computation cost.
Example: Convex ERM.
For convex ERM, we can use [15]} to instantiate the oracle. In this case, accuracy \(\alpha _{\text{private}}\)  is the excepted excess empirical risk, which is \(\alpha _{\text{private}}(n,\epsilon ,\delta ) = O{\frac{GD\sqrt{d}\sqrt{\log {1/\delta }}}{n \epsilon }}\) . Using alg:interpolate, given \(0\le \rho \le 1\) , at any point in the stream, we have,
\({\hat{F}_S(\hat{}_S)- \hat{F}_S(^*_S)} &\le \alpha _\text{private}{\frac{n}{2},\rho \epsilon ,\delta _{\text{group}}{\left\lfloor \frac{1}{\rho }\right\rfloor , \rho \epsilon , \rho \delta }} \\&\le O{\frac{GD \sqrt{d} \sqrt{\log {(1/\rho \delta ) \exp {\rho \epsilon {\left\lfloor \frac{1}{\rho }\right\rfloor -1}}}}}{n \rho \epsilon }} \\& \le O{GD {\frac{\sqrt{d}\sqrt{\log {1/\rho \delta }}}{n \rho \epsilon } + \frac{\sqrt{d}}{ n \rho \sqrt{\epsilon }}}}\\& \le O{\frac{GD\sqrt{d}\sqrt{\log {1/\rho \delta }}}{n \rho \epsilon }}\) 

where the last inequality holds when \(\frac{\epsilon }{\log {1/\rho \delta }} \le O(1)\) , which usually is the case in DP, and so is a reasonable regime. We now compare against [26]} - we ignore \(G,D\)  and \(\text{log}\)  factor in both the bounds. To have the same runtime, we need \(\rho k Tm = k^2 n \iff \rho = \frac{kn}{Tm} = \frac{kd}{\epsilon ^2n}\) , where in the last equality we substituted \(Tm = \frac{(\epsilon n)^2}{d}\) , parameters for the DP convex ERM algorithm. Our accuracy bound is \(O{\frac{\sqrt{d}}{n \rho \epsilon }} = O{\frac{\epsilon }{k \sqrt{d}}}\) , which is smaller than that of [26]}, when \(\frac{\epsilon }{k \sqrt{d}} \le {\frac{\sqrt{d}}{n k \epsilon }}^{2/5} \iff \epsilon ^7 \le \frac{\sqrt{d}^7 k^3}{n^2} \iff \epsilon \le \frac{\sqrt{d}k^{3/7}}{n^{2/7}}\) . Hence in regimes where the unlearning parameter \(\epsilon \)  is small enough, which corresponds to a stronger unlearning criterion, this algorithm is better than that of [26]}.

Experiments
We run experiments on MNIST [29]}, a standard digit classification computer vision dataset with 10 classes. We train a logistic regression model, which can be formulated as a smooth convex risk minimization problem.
Starting with a training dataset of 60k points, we simulate a stream of 300 deletions of randomly chosen points and 300 insertions of new points, randomly permuted.
We use Algorithm REF  as the learning algorithm, and the corresponding alg:unlearn-noisy-m-sgd as the unlearning algorithm.
We train for \(T=200\)  iterations, with mini-batch of size \(m=50\)  with a constant learning rate \(\eta = 0.05\) .
We run experiments on a range of values of standard deviation \(\sigma \)  of Gaussian noise, from 0 to \(1.1\)  separated by the intervals of size \(0.005\) .
For every value of \(\sigma \) , we run 10 instances of the whole unlearning procedure and report average performance: accuracy and number of unstable edits (i.e. number of times a recompute is triggered) , and their standard deviations.
Note that \(\sigma =0\)  corresponds to standard mini-batch SGD, and therefore the accuracy obtained is the accuracy for the standard method with the aforementioned setting of the hyperparameters. Moreover, the \(\sigma =0\)  setting also corresponds to Algorithm REF , and therefore the corresponding unlearning algorithm alg:unlearn-subsample-gd handles edits for this case.
In fig:exp1, we report the test accuracy (fraction of mis-classified samples in the test set) and the number of unstable edits i.e the number of times a retrain is triggered, as a function of \(\sigma \) . As expected, as \(\sigma \)  increases, we get less unstable edits. Interestingly, for small values of \(\sigma \) , like \(0.1\) , the degradation in accuracy is not as much as compared to decrease in the the number of unstable edits. Furthermore, recall that the unlearning algorithm triggers a partial recompute - fig:exp2 plots the average number of iterations done after an unstable edit compared to the number of iterations for a full recompute.
<FIGURE>