Beyond gradient-based ERM:
We also consider an approximate notion of unlearning, based on differential privacy, which has appeared in the literature [1]}, [2]}. With such a notion, we show a simple reduction to a DP algorithm, to handle unlearning requests, and show how to use group privacy to trade-off accuracy and runtime. For convex ERM, this method performs competitively with existing works (see sec:approx-unlearning).
Our proofs are simple and conceptual and so we present the key ideas in the main text and defer the proofs to the appendix.
Related work
The problem of exact unlearning in smooth convex ERM has not been studied before, and therefore the only baseline is re-computation (using some variant of gradient descent).
The most related are the works of [3]} and [1]}, which we discuss as follows.
[3]} studied the problem of \(k\) -means clustering with exact unlearning in a streaming setting of deletion requests - we borrow the setting (while also allowing insertions) and the notion of exact unlearning from therein.
We note that in [3]}, the notion of efficiency is based on the amortized (over edits) unlearning time being at most the training time since that is a natural lower bound on the overall computational cost. We, on the other hand, do not place such a restriction and so our methods can have unlearning runtime smaller than the training time. Most importantly, the general framework here (of TV-stable methods and coupling based unlearning) captures the quantized-\(k\) -means algorithm of [3]} as a special case (see sec:qsgd for details).
The work of [1]} focuses on unlearning in convex ERM problems, with a stream of edit requests, the same as here.
However there are two key differences. First, the notion of unlearning in [1]} is approximate, based on \((\epsilon , \delta )\) -differential privacy, whereas we focus on exact unlearning. Second, the unlearning runtime in [1]} is deterministic, whereas that of ours is random. These are akin to Monte-Carlo vs. Las Vegas-style of guarantee discrepancy.
We refer the reader to an extended literature survey along with a detailed comparison to [1]} in sec:detailed-related-work.
We show therein that with the same unlearning time, the accuracy guarantees of [1]} are better than us only in regimes where their approximate unlearning parameters and hence the notion, is very weak.

Problem setup and preliminaries
Streaming edit requests and exact unlearning
We describe the setup very generally.
Let \(\)  be the data space, \(\Theta \)  the output/parameter space, and \(\)  be the metadata/state space, which will be made clear later.
A procedure is a tuple \((\mathbf {A}(\cdot ),\mathbf {U}(\cdot ))\) , where \(\mathbf {A}: ^* \rightarrow \Theta \times \)  is the batch learning algorithm, and \(\mathbf {U}:\Theta \times \times \rightarrow \Theta \times \)  is the unlearning algorithm which updates the current model (first argument) and meta data (second argument) given an edit request (third argument).
Examples of meta-data could be a compressed sketch of the data points, or intermediate computations/state, which could be used upon edit time.
Let \((\cdot )\)  denote the first output of \(\mathbf {A}\)  i.e. \((\cdot ) = \mathbf {A}_1(\cdot )\) . Similarly, let \((\cdot )\)  denote the first output of \(\mathbf {U}\) .
We remark that when we refer to the algorithm's output, we mean the model output and does not include the metadata.
Finally, given two sets \(S\)  and \(S^{\prime }\) , we define \(\Delta (S,S^{\prime })\)  to be the symmetric difference between these sets i.e. \(\Delta (S,S^{\prime }) = {S \backslash S^{\prime }}+{S^{\prime } \backslash S}\) . We now define exact unlearning.
[Exact unlearning]

We say a procedure \((\mathbf {A},\mathbf {U})\)  satisfies exact unlearning  if for any
\(S, S^{\prime } \subset ^*\)  such that \(\Delta (S,S^{\prime })=1\) ,
\(\mathbf {A}(S^{\prime }) = \mathbf {U}(\mathbf {A}(S),S^{\prime }\backslash S \cup S\backslash S^{\prime })\) .
For randomized procedures, we want that for any measurable event \(\mathcal {E} \subseteq \Theta \times \) , we have \(¶{\mathbf {A}(S^{\prime }) \in \mathcal {E}} = ¶{\mathbf {U}(\mathbf {A}(S),S^{\prime }\backslash S \cup S\backslash S^{\prime }) \in \mathcal {E}}\)

A relaxation of the above definition is to maintain that only the output and not the meta-data satisfy the above condition i.e. \((S^{\prime }) = (\mathbf {A}(S),S^{\prime }\backslash S \cup S\backslash S^{\prime },S)\) . However, we will work with the stronger notion. This, with a slight difference, is referred to as perfect unlearning in [1]}.

Even though the above definition is for one edit request, it can be generalized for a stream of \(k\)  edit requests, by having that this condition holds inductively for every point in the stream.

Let \(S=S^0={_1,_2,\ldots , _n}, _i \in \)  be the initial dataset. We observe \(k\)  edit requests, each being either an insertion or deletion request. We use \(S^i\)  to denote the set of data points available at time \(i\)  in the stream. For notational simplicity, as in [1]}, we assume that at any point in the stream, the number of available data points is at least \(n/2\)  and at most \(2n\) .

Convex risk minimization
We recall some basics from convex optimization.
Let \(\subset ^d\)  be a closed convex set such that \(\text{diameter}()\)  \(\le D\)  where the diameter is measured in Euclidean distance.
Let \(\)  be the instance space and let \(f:\times \rightarrow \)  be an \(L\) -Lipschitz convex function in its first argument.
For the constraint set \(\) , given a point \(\) , a projection function \(:^d \rightarrow ^d\)  returns \(() \in \arg \min _{}{-.The function f is L-smooth in its first argument if {\nabla _f(_1,) - \nabla _f(_2,)}\le L{_1-_2} \ \forall _1,_2 \in , \in .We will drop the subscript  in \nabla  from here on.In this work, we will be concerned with smooth Lipschitz convex functions.}\)
Empirical Risk Minimization (ERM).
Given data points \(S={_1,_2,\ldots , _n}\) , we look at the following class of problems, known as empirical risk minimization (ERM).
\(\min _{\in }{\hat{F}_S():=\frac{1}{n} \sum _{j=1}^n f(,_i)}\) 
Let \((S)\)  be the output of algorithm \(\)  on dataset \(S\) . We will give guarantees on expected excess empirical risk, which is \({\hat{F} ((S))} - F(_S^*)\) , where \(^*_S\)  is the minimizer: \(^*_S \in \arg \min _{\in } \hat{F}_S()\)  and the expectation is taken with respect to the randomness in algorithm \(\) . Finally, we note that the above class of problems is large enough to be applicable even beyond machine learning - for example: in statistics, estimation problems which reduce to optimization problems of the above form are called \(M\) -estimation.
We discuss the related notion of population risk and the problem of risk minimization along with our results in sec:pop-risk.

Total variation stability and maximal coupling
We first state the definition of Total Variation (TV) distance between two distributions \(P\)  and \(Q\) .
\(\text{TV}(P,Q) = \sup _{\text{measurable sets }R} {P(R)-Q(R)} = \frac{1}{2}{\phi _P - \phi _Q}_1\) 

where the second equality holds if both distributions have probability densities with respect to a base measure which are denoted by \(\phi _P\)  an \(\phi _Q\)  respectively.
We now define total variation stability (TV-stability), which is the notion of algorithmic stability we will use.
[\(\rho \) -TV-stability]
An algorithm \(\)  is said to be \(\rho \) -TV-stable if
\(\sup _{S, S^{\prime }:\Delta (S,S^{\prime })= 1} \text{TV}((S),(S^{\prime }))\le \rho \) 

The above definition of TV stability considers the marginals of output and does not include the metadata.

Suppose \(S\)  is a dataset of \(n\)  points, and \(S^{\prime }\)  is a dataset of \(n+k_2\)  points such that \({S \backslash S^{\prime }} = k_1\) . Then, if algorithm \(\)  is \(\rho \) -TV stable, then by triangle inequality of TV and repeated applications of the above definition, we have that \(\text{TV}((S),(S^{\prime }))\le (2k_1+k_2)\rho \)

We discuss the maximal coupling characterization of total variation distance, which is a key ingredient in the design of our unlearning algorithms.
Coupling and total variation distance:
A coupling between two probability distributions \(P\)  and \(Q\)  over a common measurable space \((,)\) , where \(\)  denotes the sigma-algebra on \(\) , is a distribution \(\pi \in (\times , \otimes )\)  such that the marginals along the projections \((x,y)\rightarrow x\)  and \((x,y)\rightarrow y\)  are \(P\)  and \(Q\)  respectively. Let \(\Pi (P,Q)\)  denotes the set of couplings between \(P\)  and \(Q\) . The following describes the maximal coupling characterization of total variation distance.

For any coupling \(\pi \in \Pi (P,Q)\) , if the random variable \((p,q)\sim \pi \) , then \(\text{TV}(P,Q)\le ¶{p\ne q}\) .

There exists a “maximal" coupling \(\pi ^*\)  such that if \((p,q)\sim \pi ^*\) , then \(\text{TV}(P,Q)= ¶{p\ne q}\)

The above establishes that \(\text{TV}(P,Q)= \inf _{\pi \in \Pi (P,Q)}\mathbb {P}_{(p,q)\sim \pi }\left[p\ne q\right]\) .
As a final remark, in this work, we routinely deal with distances and divergence between probability distributions. In some cases, we abuse notation and write a divergence between random variables instead of probability distributions - these should be interpreted as the law of the random variables.

Main results
We state our main result on designing learning and unlearning algorithms in a stream of edit requests.
[Main Theorem]
For any \(\frac{1}{n} \le \rho < \infty \) , there exist a learning and a corresponding unlearning algorithm such that for any \(f(\cdot , )\) , which is \(L\) -smooth and \(G\) -Lipschitz convex function \(\forall \ \) , and a stream of edit requests,

Satisfies exact unlearning at every time point in the stream of edit requests.

At time \(i\)  in the stream, outputs \(\hat{}_{S^i}\)  with excess empirical risk bounded as,
\({\hat{F}_{S^i}(\hat{}_{S^i}) - \hat{F}_{S^i}(_{S^i}^*)} \lesssim \min {\frac{GD}{\sqrt{\rho n}},{\frac{L^{1/4}GD^{3/2} \sqrt{d}}{(\rho n)}}^{4/5}}\) 

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

We make some remarks about the result.
Training time:
Informally, what the above theorem says is that the algorithms satisfy exact unlearning and are accurate while only recomputing a \(\rho \)  fraction of times - this is indeed the nature of our algorithms.
Therefore, "Training time" here refers to the runtime of the learning algorithm.
If we measure training time in terms of number of gradient (oracle) computations, as is typical in convex optimization, then for the above accuracy, our algorithm has optimal oracle complexity in most regimes (see details in sec:learning-runtime).

Role of \(\rho \) :
The external parameter \(\rho \)  controls the trade-off between accuracy and unlearning efficiency. In the extreme case where we don't care about unlearning efficiency and are fine with paying retraining computation for every edit request, then we can set \(\rho >1\)  as large as we want to get, as expected, arbitrary small excess empirical risk. However, the interesting case is when we set \(\rho <1\) : herein, we get an improved (see below) unlearning time and yet a non-trivial accuracy, upto \(\rho \gtrsim \frac{1}{n}\) .

Strict improvement:
The above result may seem like a trade-off, but, as we argue below, is a strict improvement over the baseline of retraining after every edit request (which is the only other known method for exact unlearning for this problem). Let the target excess empirical risk be \(\alpha > \alpha _0 =\min {\frac{GD}{\sqrt{ n}},{\frac{L^{1/4}GD^{3/2} \sqrt{d}}{n}}^{4/5}}\) . For any such \(\alpha \) , there exists a \(\rho <1\) , such that our algorithms have \(\rho k \cdot \text{Training time}(\alpha )\)  expected unlearning time, which is smaller than \(k\cdot \text{Training time}(\alpha )\)  - the cost of retraining after every edit request.
Furthermore, as remarked above, since our training time is optimal in number of gradient computations (for the said accuracy), the aforementioned improvement holds for re-computation with any first-order optimization algorithm.
A small caveat is that we are comparing our expected unlearning time with deterministic runtime of retraining. To summarize, with this caveat, we have a strict improvement in the low accuracy regime, whereas in the high accuracy regime: \(\alpha < \alpha _0\) , our unlearning algorithms are as good as trivial re-computation.
However, this low accuracy regime is often the target in machine learning.
To elaborate, the goal is to minimize the population risk rather than empirical risk, and it is well known that this statistical nature of the problem results in an information-theoretic lower bound of \(\frac{1}{\sqrt{n}}\)  on excess population risk. We show in sec:pop-risk that our algorithm guarantees an excess population risk of \(\frac{1}{\sqrt{n}}+\alpha \) , and so a very small \(\alpha \)  only becomes a lower order term in excess population risk.

Algorithms:
The first upper bound on accuracy in thm:main-result is obtained by standard SGD, which, in each iteration samples a fraction of datapoints, called mini-batch, to compute the gradient, and performs the descent step - we call this sub-sample-GD.
The second upper bound is obtained using noisy accelerated mini-batch-SGD (noisy-m-A-SGD), which is also used for differentially private ERM.
Our unlearning algorithm for sub-sample-GD is rather straightforward, and most of the work is design of unlearning algorithm for noisy-m-A-SGD, which is based of efficient coupling of Markov chains corresponding to the learning algorithm. We describe the algorithms in detail in sec:algorithms.

Sub-optimality within the TV stability framework:
If we consider \(L,G,D=O(1)\) , and a simple model of computation wherein we pay a unit computation when we recompute, otherwise not, then the unlearning problem is equivalent to design of \(TV\) -stable algorithms, and a corresponding (maximal) coupling (see sec:optimal-transport for more details). Our coupling construction for unlearning in noisy-m-A-SGD, though efficient, is not maximal - this gap shows up in the accuracy bound (second term), which is \({\frac{\sqrt{d}}{\rho n }}^{1/5}\)  worse than what we would have obtained via a maximal coupling i.e \(\frac{\sqrt{d}}{\rho n}\) .
We also note that in case we don't use acceleration, but rather vanilla noisy mini-batch SGD, and the "same" coupling construction for unlearning, then we obtain a worse accuracy bound of \({\frac{\sqrt{d}}{\rho n }}^{2/3}\)  (see sec:noisy-m-sgd for details).
Finally, apart from closing the gap with the maximal coupling, another potential improvement is by giving \(\rho \) -TV stable algorithms with better accuracy. We discuss such upper and lower bounds as follows.
As pointed out, intermediate to the result in thm:main-result is the design and analysis of \(TV\) -stable algorithms for smooth convex ERM. Our main result on upper bounds on accuracy of such algorithms is the following.
[Upper bound]

For any \(0<\rho <\infty \) , there exists an algorithm which is \(\min {\rho ,1}\) -TV stable, such that for any \(f(\cdot , )\)  which is \(L\) -smooth and \(G\) -Lipschitz convex function \(\forall \ \) , and any dataset \(S\)  of \(n\)  points, outputs \(\hat{}_S\)  which satisfies the following.
\({\hat{F}_S(\hat{}_S) - \hat{F}_S(^*_S)} \lesssim GD\min {\frac{1}{\sqrt{\rho n}}, \frac{\sqrt{d}}{\rho n}}\) 
We show that the condition \(\rho \ge \frac{1}{n}\)  in thm:main-upper-bound is fundamental for any non-trivial accuracy, as evidenced by our lower bounds, with a matching dependence on \(\rho \) . Furthermore, we omit the regime \(\rho \ge 1\)  in our lower bound since it puts no constraint on the algorithm.
[Lower bound]

For any \(\rho \) -TV-stable algorithm \(\) , there exists a \(G\) -Lipschitz convex function \(f\)  and a dataset \(S\)  of \(n\)  points such the expected excess empirical risk is lower bounded as:

For any \(0<\rho <1\) , and any dimension \(d\) , \({\hat{F}_S((S)) - \hat{F}_{S}(_{S}^*)} \gtrsim GD\min {1,\frac{1}{\rho n}}\)

Assuming that \((S)\)  has a probability density function upper bounded by \(K \le O(2^d)\) , then for \(n>72,\frac{1}{n}\le \rho \le \frac{1}{4}\)  and large enough \(d\) , \({\hat{F}_S((S)) - \hat{F}_{S}(_{S}^*)} \gtrsim GD\min {{1,\frac{1}{\sqrt{\rho n}}}}\)

In each of the lower bounds, the term \(GD\)  is trivial as it is attained if an algorithm outputs a constant regardless of the problem instance.
The first lower bound holds for all problem instances without any assumptions on the relationship between the problem parameters \(d\) , \(n\)  and \(\rho \) .
Note that if we assume that the upper bound given by thm:main-upper-bound were tight, in that case we would expect to derive a lower bound of
\(\frac{\sqrt{d}}{\rho n}\)  whenever \(\frac{\sqrt{d}}{\rho n} \le \frac{1}{\sqrt{\rho n}} \iff d \le \frac{1}{\rho n}\)  - we would therefore need to shrink the class of problem instances explicitly.
Unfortunately, our techniques currently do not show improvement with this restriction.
The second result is obtained by a direct analysis, where the key ingredient is the fact that normalized volume of spherical cap of a hypersphere goes to 0 as \(d\rightarrow \infty \) , for a fixed width of the cap.
The condition that the probability distribution \((S)\)  has bounded density prevents it to have discrete atoms - this is not desirable especially since our (upper bound) algorithm sub-sample-SGD outputs a mixture of discrete distributions, and therefore does not lie in this class.
Please see sec:lower-bounds for derivations of the lower bounds.

Main ideas
In this section, we discuss the key ideas to our approach. The first is identifying a notion of stability. For this, we connect the problem of unlearning to optimal transport, and specifying a simple model of computation, the notion of total variation stability arises naturally. The second is the design of TV stable algorithms for convex ERM.
One of the algorithms we propose is an existing differential private solution, which we show to be TV stable as well. We also discuss certain important differences between our setup and that of differential privacy. Finally, the bulk of the work, is the design and analysis of efficient unlearning algorithms. We show that this problem can be reduced to efficiently constructing couplings between Markov chains, and we give such a construction using rejection sampling and reflection mappings. We now discuss these one by one.
Total variation stability from optimal transport
In this section, we give a didactic treatment of our approach to motivate the notion of total variation stability.
Consider neighbouring datasets \(S\)  and \(S^{\prime }\)  and let \(P = (S)\)  and \(Q=(S^{\prime })\)  for some randomized algorithm \(\) . The algorithm first computes on \(S\) , and then observes edit requests which generate \(S^{\prime }\)  as the current dataset. To satisfy exact unlearning , we need a procedure which moves \(P\)  to \(Q\) .
This is akin to the well-studied optimal transport problem [15]}, which briefly explain below.
Given probability distributions \(P\)  and \(Q\)  over measurable space \(\) , and a cost function \(c: \times \rightarrow \) , the goal is to transport from \(P\)  to \(Q\)  using the minimum cost.
Formally, let \(\Pi (P, Q)\)  denote the set of couplings (or transport plans) of \(P\)  and \(Q\) ; the modern (Kantorovich's) formulation asks for a transport plan \(\pi \)  which minimizes the expected cost: \( \min _{\pi \in \Pi (P,Q)} \mathbb {E}_{(x,y)\sim \pi } c(x,y)\) .
A model of computation:
Note that there is of course the trivial coupling in which we generate independent samples from \(P\)  and \(Q\)  - this corresponds to re-computation, which, as argued, is not an efficient method in general.
Instead, we should correlate \(P\)  and \(Q\)  so that transporting from \(P\)  to \(Q\)  can reuse the randomness (computation) used for \(P\) .
For this, we use the cost function in the optimal transport problem as a surrogate of modelling computation.
In the optimal transport problem, the cost is typically a distance on the space, whereas we are concerned with computational cost. So is there a distance function which corresponds to computational cost? Note that
the sequential nature of the problem already gives us samples generated from \(P\) , so a natural question is, can we use this to transport to \(Q\) ?
We can set the cost function as \(c(x,y)= {\left\lbrace \begin{array}{ll} 1 & \text{ if } x\ne y \\ 0 & \text{ otherwise}\end{array}\right.}\) .
This corresponds to an oracle which charges a unit computation if we use \(y\)  which is different from \(x\) , which can correspond to a recomputation. Under this simple model of computation, the optimal expected computational cost becomes exactly equal to the total variation distance between \(P\)  and \(Q\) :
\(\inf _{\pi \in \Pi (P,Q)} \mathbb {1}{x \ne y}\)  - the maximal coupling characterization of total variation distance.

TV stability:
The above establishes that if we want to transport \(P\)  to \(Q\)  using minimum computation cost, the expected computation cost cannot be smaller than the total variation distance between \(P\)  and \(Q\) .
Intuitively, this means that is least \(1-\text{TV}(P,Q)\)  fraction of samples are representative for both \(P\)  and \(Q\) .
From the sequential nature of our problem, when we generate \(P\)  - the output on dataset \(S\) , we don't know what \(Q\)  would be, since we don't know the incoming edit request. Hence a reasonable property to have in the algorithm is that its output is close in total variation distance uniformly over all possible \(Q\) 's. This motivates our definition of total variation stability.

Optimal transport vs unlearning:
Unlike the optimal transport problem wherein we are given \(P\)  and \(Q\) , and the task is to find a coupling, in our setup, we have to find an algorithm generating \(P\)  and \(Q\)  as well as the coupling.
Moreover, for a fixed \(\rho \) , there may be many algorithms which are \(\rho \) -TV stable.
The goal therefore, is to find among these algorithms, the one with the maximum accuracy for the (convex ERM) problem, and for which we can design a corresponding efficient unlearning algorithm.

TV-stable learning algorithms and differential privacy
In this section, we discuss the ideas underlying the design of TV-stable learning algorithms. We first give the definition of differential privacy (DP), which will be a key tool.
Differential privacy is a notion of data privacy, introduced in [16]}, defined as follows.
[Differential privacy (DP)]

An algorithm \(\mathcal {A}\)  satisfies \((\epsilon ,\delta )\) -differential privacy if for any two neighboring datasets \(S\)  and \(S^{\prime }\) , differing in one sample, for any measurable event \(\in \text{Range}()\) ,
\(\mathbb {P}(\mathcal {A}(S)\in ) \le e^\epsilon \mathbb {P}(\mathcal {A}(S^{\prime })\in )+\delta \) 


Intuitively, a differentially private algorithm promises that the output distributions are close in a specific sense: the likelihood ratios for all events for two neighbouring datasets is uniformly close to \(e^{\pm \epsilon }\) , upto a failure probability \(\delta \) .
Note that we have identified that we want our outputs to be \(\rho \) -TV stable. A natural question is whether we can relate the \((\epsilon , \delta )\) -DP notion to \(\rho \) -TV-stability. An easy to see direction is that any \(\rho \) -TV stable method is (at least) \((0,\rho )\) -DP.
Similarly, for the other direction, under additional assumptions, such relations can be derived. The important part is that certain widely used DP methods are TV stable as well. The primary example, which we will use in this work, is Gaussian mechanism. It is known that adding Gaussian noise of variance \(\frac{\sqrt{\log {1/\delta }}}{\epsilon }\)  to a 1-sensitive function, provides \((\epsilon ,\delta )\) -DP [17]}. It can be shown that the same method also provides \(\rho \) -TV stability, with \(\rho = \frac{\epsilon }{\sqrt{\log {1/\delta }}}\) .
TV-stable algorithm:
For the problem of TV stable convex empirical risk minimization, we propose two algorithms: sub-sample-GD and noisy-m-A-SGD. We show that the expected excess empirical risk of noisy-m-A-SGD is better than that of sub-sample-GD, in regimes of small dimension.
The algorithm noisy-m-A-SGD is essentially noisy-m-SGD algorithm, which appeared in [18]} for DP convex ERM, with an additional Nesterov's acceleration on top.
For ease of presentation, we discuss using noisy-m-SGD (no acceleration), and in
paragraph titled "Fast algorithms and maximal coupling" in the next section, it would become clear why adding acceleration helps us.
In noisy-m-SGD,
at iteration \(j\) , we sample a mini-batch \(b_j\)  uniformly randomly, use it to compute the gradient on the previous iterate \(_{j}\)  denoted as \(\nabla \hat{F}_S(_{j},_{b_{j}})\)  and update as follows:
\(_{j+1} = _j - \eta {\nabla \hat{F}_S(_{j},_{b_{j}}) + \theta _t}\) 
where \(\theta _j \sim (0,\sigma ^2_d)\)  and \(\sigma \)  is set appropriately.
We ignore the projection step in the current discussion.
This procedure can be viewed as sampling from a Markov chain depicted in Figure REF .
<FIGURE>
DP convex ERM and unlearning:
We first discuss an important distinction in the differential privacy and our unlearning setup.
In the DP setup, we have a curator which possesses the dataset, and an analyst/adversary, against which the curator want to provide privacy. The analyst queries the dataset, and the curator provides DP answers to the queries. The curator can also reveal additional information pertaining to the algorithmic details, however, it is beneficial to the curator to only release limited information. In particular, the curator can chose to keep certain states of the algorithm secret.
This could be done in the case when only the marginals of the output satisfy a strong DP-guarantee.
So, if the curator were to release the secret state as well, the adversary can correlate information and then the privacy level, which is now measured using the joint distribution of output and state, degrades.
In the noisy-m-SGD algorithm for example, the output typically is the average or final iterate whereas the rest of iterates and mini-batch indices \(b_j\) 's are the secret state.
In the unlearning setup, there is no adversary per se, or in the idealized application, the curator is the adversary and the dataset owners wants it to have as little control as possible.
It is therefore natural to
demand that the probability distribution of the entire state maintained by the algorithm, and not just the output be exactly identical after performing the unlearning operation.
This, with slight differences, is referred to as perfect unlearning in [1]}, and what our algorithms satisfy.
We have argued that designing TV stable algorithms is a good start, and for a moment suppose that the TV stability is same as DP.
Then should we measure TV stability between the joint distributions over the entire state?
This would limit the application of DP techniques in which keeping additional state hidden has stronger privacy property.
In that case, TV stability parameter, and hence the computational cost of unlearning would be large.
Interestingly, even though the previous work in differentially private convex ERM, for example [18]}, argues that the released iterate (average/final iterate) is differentially private, the analysis is typically carried out by first arguing, via a composition step, that all iterates together are differentially private. This means that all iterates can be released without any additional cost of privacy.
This innocuous property arguably provides no benefit for privacy, but turns out to be extremely beneficial to us in unlearning. However, even though the all the iterates can be released, the mini-batches still need to kept secret.
We handle this in the unlearning algorithm using an estimation step - see paragraph titled “Estimation of marginals" in sec:main-idea-unlearning.

Unlearning via (un)couplings
The final, though the most important piece, is the design of unlearning algorithms.
Recall that \(S\)  is the initial dataset, \(S^{\prime }\)  is the dataset after one edit request, and we want to design a transport from \(P=(S)\)  to \(Q=(S^{\prime })\) , which means that we need to construct a coupling of \(P\)  and \(Q\) .
Broadly, there are two challenges: the first is the data access restriction - when generating a sample from \(P\) , we don't know what \(Q\)  would be, therefore, the coupling cannot be based on efficiently sampling from a joint distribution directly, but is limited to work with samples generated from \(P\) .
The other is that construction of the coupling should be computationally more efficient than drawing independent samples from \(P\)  and \(Q\) , which essentially amounts to our baseline of re-computation.
An efficient general approach:
We first setup some terminology - the diagonal of a coupling \(\pi \)  of two probability distributions, is the set \({(p,q): p = q}\)  where \((p,q)\sim \pi \) , and similarly, the non-diagonal is the set \({(p,q): p \ne q}\) , \((p,q)\sim \pi \) . We have that the measure of the non-diagonal, under a maximal coupling \(\pi ^*\) , is \(\mathbb {P}_{(p,q)\sim \pi ^*} \mathbb {1}{(p,q): p \ne q} = \text{TV}(P,Q)\) .
This implies that when using \(\rho \) -TV stable algorithms, the probability measure of the diagonal under a maximal coupling, is large - at least \(1-\rho \) .
At a high-level, our unlearning approach comprises of two stages: verification and recomputation.
We first verify whether our output on dataset \(S\)  (i.e. sample from \(P\) ) falls on the diagonal of any maximal coupling of \(P\)  and \(Q\)  or not - if that is indeed the case, then the same sample for \(Q\)  suffices.
For computational efficiency, we require that verification be computationally much cheaper then recomputing (smaller that \(\rho \cdot \text{recompute cost}\) ).
If the verification fails, we sample from the non-diagonal of any maximal coupling \(P\)  and \(Q\) , so that we have a valid transport.
As the name suggest, the computational cost of recomputation that we will shoot for is to be of the same order as (full) recompute.
If we are able to design such a method, then we will show that for \(k\)  edit requests, the expected computational cost for unlearning is \(k \cdot \text{verification cost} + k\rho \cdot \text{recompute cost} \approx k\rho \cdot \text{recompute cost}\) .
The design of unlearning algorithm is dependent on the corresponding learning algorithm.
Since we proposed two learning algorithms, each has a corresponding unlearning algorithm. Herein, we will only discuss the more challenging case, which is for the noisy-m-SGD (no acceleration) algorithm.

Coupling of Markov chains:
Our approach for unlearning is to construct a coupling of the optimization trajectories on neighbouring datasets.
We have discussed that the iterates from noisy-m-SGD can be seen as generated from a Markov chain, depicted in Figure REF .
Hence, for two neighbouring datasets, the iterates are sampled from two different Markov chains \(P\)  and \(Q\) .
Moreover, by design, we know that these Markov chains are \(\rho \) -TV close - we measure the total variation distance between joint distribution of marginals of iterates i.e. \(\text{TV}{{_j^P}_{j=1}^{T+1},{_j^Q}_{j=1}^{T+1}} \le \rho \) .
The task is now to maximally couple these two Markov Chains.
We remark that in the Markov chain literature, maximal coupling of Markov chains does not refer to the above but rather the setting wherein we have one Markov chain, but started at two different states, and the goal is to design a coupling such that their sampled states become and remain equal as soon as possible. In contrast, our notion of coupling of two Markov chains has also been recently studied by [21]} and [22]}, wherein they refer to this problem as design of uncoupling or maximal agreement/exit couplings.

Overview of unlearning algorithm:
The learning algorithm saves all states depicted in Figure REF . In the unlearning algorithm, we proceed sequentially: in iteration \(j\) , we first couple the mini-batches which amounts to replacing the deleted point by a uniformly random point, or inserting the new point in some mini-batches: let the coupled mini-batches be \({b^P_j}_{j=1}^T\)  and \({b^Q_j}_{j=1}^T\) .
We then compute an estimate of marginal densities of \(_j\)  under \(P\)  and \(Q\) , via conditional densities under coupled mini-batches i.e. we compute \(\phi _P(_j|b_j^Q)\)  and \(\phi _Q(_j|b_j^P)\)  - note that these are just Gaussian densities evaluated at \(_j,b_j^P\)  and \(_j,b_j^Q\)  respectively, where \(_j\)  is sample from \(P\) .
We then do a rejection sampling step wherin we draw a uniform random variable \(u \sim \text{Unif}(0,1)\) , and then check if \(u \le \frac{\phi _Q(_j|b_j^Q)}{\phi _P(_j|b_j^P)}\) .
If the step succeeds, we accept \(_j\)  as a sample from \(Q\)  and move to the next iteration and repeat. If any of the rejection sampling step fails, say at step \(t\) , we generate \(_{t+1}^Q\)  by reflecting \(_{t+1}^P\)  about the mid-point of means of the two Gaussians at step \(t\)  for \(P\)  and \(Q\) . After this reflection step, we abandon the rest of iterates from \(P\)  and generate the new iterates from \(Q\)  by continue retraining on dataset \(S^{\prime }\) .
This procedure is described as alg:unlearn-noisy-m-a-sgd (please see sec:unlearn-noisy-m-a-sgd for more details).
In the above, the rejection sampling steps comprise the verification stage, and if any of the rejection sampling fails, we move to re-computation.
The reason why verification can be done efficiently here is due to the finite sum structure of the ERM problem. To elaborate, at any iteration,
to compute the conditional density \(\phi _Q(_j\vert b_j^Q)\) , we need to compute the gradient with the new dataset \(S^{\prime }\)  - this, using the gradient of the old dataset only requires subtracting the gradient at the deleted point, so \(O_d(1)\)  runtime as opposed to \(O_d(m)\) , if we were to compute from scratch, where \(m\)  is the mini-batch size. Moreover, throughout verification, this computation is done only for iterations which used the deleted point which are roughly \(\frac{Tm}{n}\)  iterations. Hence the total runtime of verification is \(O_d\left(\frac{Tm}{n}\right)\)  as opposed to \(O_d\left(Tm\right)\)  for re-computation.
Finally, if the probability of recompute is \(\rho \) , then expected unlearning time for \(k\)  edits is \(\approx \rho kTm + \frac{kTm}{n}\)  - note that the second (verification time) is a lower-order term as long as \(\rho \gtrsim \frac{1}{n}\)  - furthermore, this \(\rho \gtrsim \frac{1}{n}\)  is the best possible within the TV stability framework for non-trivial accuracy as evidenced by our lower bounds (see thm:main-lower-bound).
Please see sec:unlearning-runtime for more details on runtime and efficient implementation using suitable data structures.

Fast algorithms and maximal coupling:
The above procedure generates a coupling but not a maximal coupling - the measure of the diagonal under the coupling, and hence the probability to recompute, is \(\sqrt{T}\)  worse then the optimal, where \(T\)  is the number of iterations run of noisy-m-SGD. This gives us that the faster the algorithm (in terms of iteration complexity) is, the smaller the probability to recompute, when using our coupling construction. This motivates why we use accelerated mini-batch SGD, since it has a quadratically faster iteration complexity than vanilla mini-batch SGD. In sec:learning-runtime, we also remark that using even faster algorithms like Katyusha [23]} does not yield further improvements. Finally, the design of maximal coupling would (likely) be done via one step rejection sampling, instead of doing it iteratively. However, if the rejection sampling fails, sampling from the non-diagonal efficiently is tricky. We leave the question of obtaining a maximal coupling for future work.

Estimation of marginals:
We remarked that we want to create maximal coupling of marginals of the output, and therefore measure TV distance between marginals, rather than the entire algorithmic state.
Consider one, say \(j^{\text{th}}\) , iteration of noisy-m-SGD,
then \(b_j\)  is additional state, and we measure TV between marginals of \(_j\)  and \(_j^{\prime }\) .
The distribution of \(_j\)  is such that, for any event \(E\)  in range of \(_j\) , \(P(_j \in E) = \mathbb {E}_b P(_j \in E | b_j = b)\) .
To construct a coupling between the marginals via rejection sampling, we need to evaluate the ratio of marginal densities: \(\frac{\phi _Q(_j)}{\phi _P(_j}\) . However, the marginal is a mixture distributions with large (exponential in \(m\)  (mini-batch size)) number of components, and therefore even evaluating the marginal density is infeasible. One solution is to just consider \(m=n\)  i.e. full gradient descent, and then there is no additional state.
However, this makes the training runtime worse, and that means that we would be using a slower learning algorithm than what we would have used if we were to simply recompute to unlearn.
Hence, to tackle this, as described in the previous paragraph, we evaluate the ratio of conditional probability densities, where the conditioning is on the coupled mini-batch indices (\(b_j^P\)  and \(b_j^Q\) ) i.e. \(\frac{\phi _Q(_j|b_j^Q)}{\phi _P(_j|b_j^P)}\) .
This corresponds to using unbiased estimates of the marginals densities.
It is easy to verify, using convexity of the pointwise supremum for instance, that \(\text{TV}((^P_j,b^P_j),(^Q_j,b^Q_j)) \ge \mathbb {E}_{(b_j^P,b_j^Q)} \text{TV}(_j^P|b_j^P, _j^Q|b_j^Q) \ge \text{TV}(_j^P,_j^Q)\) .
However, in general, this might still not be ideal since we are estimating with just one sample from the mixture and hence the estimation error would be large. However, we will show that since we are anyway not able to construct maximal couplings, we don't pay extra with this coarse estimate.

Algorithms
In this section, we present the algorithms for learning and unlearning. In our algorithms, we use functions “save" and “load", which vaguely means saving and loading the variables to and from memory respectively. In sec:runtime, we explain what data structures to use for computational efficiency.
The proofs of results in this section are deferred to sec:proofs-learning.
TV-stable learning algorithms
sub-sample-GD
The first algorithm, which is superior in high dimensions, called sub-sample-GD, is just vanilla mini-batch SGD wherein at each iteration, a mini-batch of size \(m\)  is sub-sampled uniformly randomly. Furthermore, we save all the mini-batch indices, gradients and iterates to memory.
We will see that the unlearning algorithm presented (alg:unlearn-subsample-gd) uses all the saved iterates. However this is done only for ease of presentation - in sec:space-complexity, we discuss a simple efficient implementation (of the unlearning algorithm), which doesn't need any iterate, yet has the same unlearning time complexity.
[h]
sub-sample-GD\((_{t_0}, t_0)\)

[1]
Initial model \(_{t_0}\) , data points \({_1,\ldots , _n}, T, m, \eta \) 
\(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)\) 
\(_{t+1} = {_t - \eta _t}\) 
Save(\(b_t,_t,_t\) )
\(\hat{}_S = \frac{1}{T}\sum _{t=1}^{T+1}_t\)
We now give guarantees on excess empirical risk for sub-sample-GD.
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\sqrt{\rho n}}{GT}}, T=\frac{DL\sqrt{\rho n}}{G}\) , and \(m= \max {{\frac{G\sqrt{\rho n}}{DL},1}}\) , outputs \(\hat{}_S\)  which is \(\min {\rho ,1}\) -TV-stable and satisfies \({\hat{F}_S(\hat{}_S) - \hat{F}_S(^*_S)} \lesssim \frac{GD}{\sqrt{\rho n}}.\)

noisy-m-A-SGD
The second algorithm, superior in low dimensions, called noisy-m-A-SGD is mini-batch accelerated SGD with appropriate Gaussian noise added at each iteration. In the literature, this algorithm (with or without acceleration) is used for DP training of (non) convex models.
In each iteration, we save the mini-batch indices, the models, the gradients as well as the noise vectors to memory.
[ht]
noisy-m-A-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
Sample \(\theta _t \sim (0,\sigma ^2_d)\) 
\(  \mathring{}_t = (1 - \alpha _t)_t + \alpha _t _{t-1} \) 
\(_t = \frac{1}{m}\sum _{j \in b_t}\nabla f( \mathring{}_t,_j)\) 
\(_{t+1} = { \mathring{}_t - \eta { _t +\theta _t}}\) 
Save(\(b_t,\theta _t,_t, \mathring{}_t, _t\) )
\(\hat{}_S =_{T+1}\)
We now state our results for noisy-m-A-SGD.
Let \(f(., )\)  be an \(L\) -smooth \(G\) -Lipschitz convex function \(\forall \ \) .
For any \(0<\rho <\infty \) ,
Algorithm REF , run with
\(t_0=1, \eta = \min {\frac{1}{2L}, \frac{D}{{\frac{G}{\sqrt{m}} +\sigma }T^{3/2}}}\) , \(\alpha _0=0, \alpha _t = \frac{1-t}{t+2}\) ,
\(\sigma = \frac{8\sqrt{T}G}{n\rho }\) , and \(T \ge \frac{(n\rho )^2}{16m^2}\)  outputs \(\hat{}_S\)  which is \(\min {\rho ,1}\) -TV stable and satisfies
\({\hat{F}_S(\hat{}_S) - \hat{F}_S(^*_S)} \lesssim \frac{LD^2}{T^2} + \frac{GD}{\sqrt{Tm}} + \frac{G D\sqrt{d}}{n \rho }.\) 
Choosing \(T\)  and \(m\)  appropriately gives us the following corollary.


Let \(f(., )\)  be an \(L\) -smooth \(G\) -Lipschitz convex function \(\forall \ \) .
For any \(0<\rho <\infty \) ,
Algorithm REF , run with \(t_0=1, m \ge \min {\frac{d}{16}, \frac{1}{4}{\frac{(\rho n)^3 G\sqrt{d}}{LD}}^{1/4}}\) ,
\(\eta = \min {\frac{1}{2L}, \frac{D}{{\frac{G}{\sqrt{m}} +\sigma }T^{3/2}}}\) , \(\alpha _0=0, \alpha _t = \frac{1-t}{t+2}\) ,
\(\sigma = \frac{8\sqrt{T}G}{n\rho }\) , and \(T=\max {\frac{(\rho n)^2}{md}, \sqrt{\frac{LD\rho n}{G\sqrt{d}}}}\)  outputs \(\hat{}_S\)  which is \(\min {\rho ,1}\) -TV stable and satisfies \({\hat{F}_S(\hat{}) - \hat{F}_S(^*_S)} \lesssim \frac{GD\sqrt{d}}{\rho n}\) .
The choice of \(T\)  in corr: upper-bound-m-A-sgd yields that the largest mini-batch size that can be set, without hurting runtime, is \(m = {\frac{(\rho n)^3 G}{\sqrt{d}^3 L D}}^{1/2} = {\frac{G}{LD}}^2T^3\) .
Furthermore, the condition \(m \ge \min {\frac{d}{16}, \frac{1}{4}{\frac{(\rho n)^3 G\sqrt{d}}{LD}}^{1/4}}\)  yields \((\rho n) \ge {\frac{LD (\sqrt{d})^7}{256 G}}^{1/3}\) .
In prop:tv-stability-noisy-sgd (in Appendix , sec:alg-learning), we show that the upper bound on total variation stability parameter of Algorithm REF 
derived in Proposition
REF  is tight in all problem parameters, upto constants.

Unlearning algorithms
We now discuss algorithms to handle edit requests which are based on efficiently constructing couplings, in some cases maximal couplings.
An important component on constructing such couplings is, what we call verification, wherein, at a high-level, we check if the current model is likely after the edit request or not. If the verification is successful, we don't do any additional computation, otherwise we do a partial or full recompute (i.e. retrain), which we call recomputation.
The key insight is that verification can be done efficiently, and fails with small probability (depending on the TV-stability parameter).
We now discuss the two algorithms, one for handling unlearning in Algorithm REF  and the other for Algorithm REF , and show that the probability with which a recompute is triggered is small - please see sec:runtime, for a finer analysis of runtime.
The proofs of results in this section are deferred to sec:proofs-unlearning.
Unlearning for sub-sample-GD
At the start of the stream, at every iteration of
sub-sample-SGD,
we sample a mini-batch of size \(m\)  out of \(n\)  points uniformly randomly, and then compute a gradient using these samples - note that this is the only source of randomness in the algorithm.
As we progress along the stream observing edit requests, the number of available data points changes.
Therefore, if the algorithm were executed on this dataset of, say \(\tilde{n}\)  points, at every iteration it would have sub-sampled \(m\)  out of \(\tilde{n}\)  (and not \(n\) ) points.
The way to account for this discrepancy is to simply adjust the sub-sampling probability measure accordingly.
Coupling mini-batch indices:
The main idea to unlearning in alg:unlearn-subsample-gd is to couple the sub-sample indices.
For deletion, we just look at each mini-batch, and (literately) verify if the deleted point were used or not. If the deletion point was not used in any iterations, then we don't do anything, otherwise, we trigger a recompute. In the case of insertion, there is no such way of selecting iterations in which the point was sampled, because the inserted point was absent. However, we know that the new point would have been sampled with probability \(m/(n+1)\) .
We can thus verify by selecting each iteration with the same probability.
We then replace a uniformly sampled point in the mini-batch of that step by the inserted point.
Algorithm REF  implements the above procedure.
[!htbp]
Unlearning for sub-sample-GD

[1]
Data point index \(j\)  to delete or data point \(\)  to insert (index \(n+1\) )
\(t=1,2\ldots , T\) 
Load\({ b_t,_t,_t}\) 
deletion and \(j \in b_t\) 
sub-sample-GD\((_t,t)\)  // Continue training on current dataset
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,)}\) 
\(_{t+1} = {_t - \eta { _t^{\prime } +\theta _t}}\) 
Save(\(_{t+1}, _t^{\prime }, b_t\backslash {i} \cup {n+1} \) )
sub-sample-GD\((_{t+1}, t+1)\)  // Continue training on current dataset
We state our main result for unlearning with alg:unlearn-subsample-gd below.
(Algorithm REF , Algorithm REF ) satisfies exact unlearning. Moreover, for \(k\)  edits, Algorithm REF  recomputes with probability at most \(2k\rho \) .

Unlearning for noisy-m-A-SGD
[ht]
Unlearning for noisy-m-A-SGD

[1]
Data point index \(j\)  to delete or data point \(\)  to insert (index \(n+1\) )
\(t=1,2\ldots , T\) 
Load\({\theta _t, _t, \mathring{}_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( \mathring{}_t,_j) - \nabla f( \mathring{}_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( \mathring{}_t,_i) - \nabla f( \mathring{}_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-A-SGD(\(_{t+1},t+1\) )  // Continue retraining on current dataset 
Our unlearning algorithm for noisy-m-A-SGD is based on efficiently constructing a coupling of Markov chain describing noisy-m-A-SGD, with large mass on its diagonal.
The key ideas have already been described in sec:main-idea-unlearning, and we just fill in some details here.
We first describe how alg:unlearn-noisy-m-a-sgd couples mini-batch indices while handling edit request.
Coupling mini-batch indices:
After observing a deletion request, in Algorithm REF ,
we look at all iterations in which the deleted point was sampled. We then replace the deleted point with a uniformly random point not already sampled in that iteration. For insertion, at each step, we again replace a uniformly sampled point in the mini-batch of that step by the inserted point with probability \(\frac{m}{n+1}\) .

Reflection maps:
We define the notion of reflection map, which will be used in our coupling construction.
[Reflection map]
Given a vector \( and two vectors \)\( and \)\(, the reflection of \)  under \((,)\) , denoted as \(\text{reflect}(,)\) , is defined as
\(\text{reflect}(,) = +(-\) 
Reflection coupling is a classical idea in probability, used to construct couplings between symmetric probability distributions [24]}.
The reflection map, given \( , \) , reflects \( about the mid-point of \)\( and \)\(.The context in which we will use it is \)  would be a sampled point from a Gaussian under old dataset \(S\)  (on which the model was trained on), and \(\)  and \(\)  being the means of the Gaussian under new dataset \(S^{\prime }\)  (after edit request) and \(S\)  respectively.
The map essentially exploits the spherical symmetry of the Gaussian to generate a good sample for the distribution under \(S^{\prime }\) . Please see sec:reflection for some properties of the reflection map, which are used in the final proofs.

Iterative rejection sampling:
Our unlearning algorithm is based on iteratively verifying each model \(_{t+1}\)  using rejection sampling. To elaborate, at each iteration, we check if the noisy iterate, defined as \(\bar{}_{t+1} =  \mathring{}_{t} - \eta (_{t} + \theta _{t})\)  is a good sample for the dataset \(S^{\prime }\) , where \(_t\)  is the gradient computed on \( \mathring{}_t\)  using a uniform sub-sample from \(S\) .
To do this, we need to compute a ratio of estimated marginal densities of \(_{t+1}\)  for both datasets, evaluated at the noisy iterate, and compare it with \(\text{Uniform}(0,1)\) .
It it succeeds, we move to the next iteration and repeat. If any of the rejection sampling fails, we do a reflection, and continue retraining on \(S^{\prime }\) .

Estimation of marginals:
We explain what we mean by estimated marginal densities in the previous paragraph.
As remarked before, if we did not sub-sample mini-batches (i.e. used gradient descent), then we would simply use the marginal distribution of iterates for rejection sampling. However, that would amount to a worse runtime. Instead, we estimate the marginal densities as follows:
fix all iterates before iteration \(t\) , and consider noisy iterate \(\bar{}_{t+1} =  \mathring{}_t - \eta (_t + \theta _t)\) .
If we also fix the sampled mini-batch \(b_t\) , then \(\bar{}_{t+1}\) 
is distributed as \(( \mathring{}_t - \eta _t, \eta ^2\sigma ^2 )\) .
However, once we unfix
\(b_t\) , then \(_{t+1}\)  is mixture of Gaussians, with the number of components being exponential in \(m\) .
Ideally, to do rejection sampling, we need to compute the marginal density of the distribution of \(_{t+1}\)  (and \(_{t+1}^{\prime }\)  - the iterate for dataset \(S^{\prime }\) ) evaluated at \( \mathring{}_t - \eta (_t+\theta _t)\)  - computing which however, is infeasible.
Therefore, we just use the coupled mini-batches indices as a sample from the mixture and estimate the marginal density using the conditional density - this is done in line 15 of alg:noisy-m-a-sgd, with a small change that we evaluate the ratio of conditional densities of noisy gradients rather than iterates, but it can be verified that the ratio is invariant to this shift and scaling.
Please see sec:unlearn-noisy-m-a-sgd-appendix for a more formal treatment of the coupling procedure. We now state the main result for this section.
(Algorithm REF , Algorithm REF ) satisfies exact unlearning. Moreover, for \(k\)  edits, Algorithm REF  recomputes with probability at most \(\frac{k\rho \sqrt{T}}{4}\)

Proofs of main results
In this section, we give the proofs of main results, stated in sec:main-results, using the results in the preceding sections.
Proof of Theorem 
The proof follows by combining the guarantees for the two algorithms we present: sub-sample-GD (Algorithm REF ) and noisy-m-A-SGD (Algorithm REF ), and their corresponding unlearning algorithms: Algorithm REF  and Algorithm REF .
We discuss these one by one.
From prop:upperboundsubsample, we have that, given \(0<\rho \le 1\) , sub-sample-GD is \(\rho \) -TV stable and has excess empirical risk bounded by \(O{\frac{GD}{\sqrt{\rho n}}}\) . This holds at every point in the stream by assumption that the number of samples are between \(\frac{n}{2}\)  and \(2n\) . Furthermore, from prop:unlearn-subsample-GD, we have that the unlearning algorithm
satisfies exact unlearning at every point in the stream, proving the first part of the claim for sub-sample-GD.
Moreover, it states that recompute probability for \(k\)  edit requests is \(O(\rho k)\) . Finally, from claim:runtimesubsampleGD, we have that there exist efficient implementations, such that the runtime of unlearning for sub-sample-GD is \(O(\max {k,\min {\rho ,1} k \cdot \text{Training time}}\) , where "Training time" is the runtime of the corresponding learning algorithm - this means that re-computations overwhelm the total unlearning time. This establishes all the guarantees for one algorithm and recovers one of the upper bounds in the second claim.
The situation for the other algorithm is a little more involved.
From prop:upper-bound-noisy-m-a-sgd, for dataset \(S\)  of \(n\)  points, we have that, given \(0<\tilde{\rho }\le 1\) , noisy-m-A-SGD is \(\tilde{\rho }\) -TV stable and its excess empirical risk is bounded as follows:
\({\hat{F}_S(\hat{}_S) - \hat{F}_S(^*_S)} \lesssim \frac{LD^2}{T^2} + \frac{GD}{\sqrt{Tm}} + \frac{G D\sqrt{d}}{n \tilde{\rho }},\) 
where \(T\)  is the number of iterations for noisy-m-A-SGD algorithm, and \(m\)  the mini-batch size.
From Proposition REF , we have that the unlearning algorithm satisfies exact unlearning (establishing the first claim) and recomputes, for \(k\)  edit requests, with probability \(O(\tilde{\rho }k \sqrt{T})\) .
Finally, from claim:runtimenoisyASGD, we have that there exist efficient implementations, such that the runtime of unlearning for noisy-m-A-SGD is \(O(\max \lbrace k,k \min {\tilde{\rho }\sqrt{T},1} \cdot \text{Training time}\rbrace )\) .
In the statement of thm:main-result, we want that the unlearning runtime be such
that we recompute for a \(\rho \)  fraction of edit requests (as opposed to something dependent on \(T\) ).
Therefore, we substitute \(\tilde{\rho }= \frac{\rho }{\sqrt{T}}\) , and this changes the excess empirical risk bound for noisy-m-A-SGD, as follows:
\({\hat{F}_S(\hat{}_S) - \hat{F}_S(^*_S)} \lesssim \frac{LD^2}{T^2} + \frac{GD}{\sqrt{Tm}} + \frac{G D\sqrt{d}\sqrt{T}}{n \rho }.\) 
We use the largest mini-batch size, which does not hurt runtime, which is \(m={\frac{G}{LD}}^2T^3\) . This simplifies the upper bound to \(\frac{LD^2}{T^2} + \frac{G D\sqrt{d}\sqrt{T}}{n \rho } \) . Optimizing the trade-off, we have \( \frac{LD^2}{T^2} = \frac{G D\sqrt{d}\sqrt{T}}{n \rho } \iff T = {\frac{LD (n\rho )}{G \sqrt{d}}}^{2/5}\) , and the excess empirical risk becomes \({\hat{F}_S(\hat{}_S) - \hat{F}_S(^*_S)} \lesssim \frac{LD^2}{T^2} = {\frac{L^{1/4}GD^{3/2} \sqrt{d}}{(\rho n)}}^{4/5}\)  – this recovers the other term in the upper bound in thm:main-result. However, note that Proposition REF  has an additional condition that \(T\ge \frac{(n\tilde{\rho })^2}{16m^2}\)  - we show that in our setting of \(\tilde{\rho }\)  and \(m\) , this condition is equivalent to the excess empirical risk of noisy-m-A-SGD being smaller than that of sub-sample-GD. Hence, the regime in which the aforementioned condition is violated is the same regime in which it is better to use the other sub-sample-GD algorithm, and therefore
is benign. Setting \(m={\frac{G}{LD}}^2T^3\)  and \(\tilde{\rho }= \rho / \sqrt{T}\) , the condition simplifies as \(T\ge \frac{(n\rho )^2}{16 T\cdot T^6 } {\frac{LD}{G}}^2 \iff T^8 \ge \frac{(n\rho )^2}{16} {\frac{LD}{G}}^4 \iff {\frac{LD (n\rho )}{G \sqrt{d}}}^{16/5} \ge \frac{(n\rho )^2}{16} {\frac{LD}{G}}^4 \iff {\frac{\sqrt{d}}{(n\rho )}}^{4/5} \le \frac{2}{\sqrt{n\rho }}{\frac{LD}{G}}^{1/5} \iff {\frac{L^{1/4}GD^{3/2} \sqrt{d}}{(\rho n)}}^{4/5} \le \frac{2GD}{\sqrt{n\rho }}\) , where the final inequality indicates that the expected excess empirical risk of noisy-m-A-SGD is at most that of sub-sample-GD, up to constants.
The above is established for dataset \(S\)  but holds for any dataset \(S^i\)  in the stream using the assumption that the number of samples are between \(\frac{n}{2}\)  and \(2n\) .
Combining the above arguments finishes the proof of thm:main-result.
\(\Box \)

Proof of Theorem  REF
We give two algorithms, sub-sample-GD (Algorithm REF ) and noisy-m-A-SGD (Algorithm REF ), one for each of the upper bounds.
From prop:upperboundsubsample and corr: upper-bound-m-A-sgd, we have that, given \(0<\rho < \infty \) , these are \(\min {\rho ,1}\) -TV stable and their excess empirical risk is bounded is \(O{\frac{GD}{\sqrt{\rho n}}}\)  and \(O{\frac{GD\sqrt{d}}{\tilde{\rho }n}}\)  respectively.
Hence combining the above by taking a minimum, establishes the claimed result.
\(\Box \)

Proof of Theorem  REF
In all the lower bounds, we have a \(GD\)  term - this is a trivial lower bound, since if an algorithm is defined as \((S)=0\)  (or any constant), then this is perfectly \(TV\)  stable (\(\rho =0\) ), and the expected excess empirical risk is upper bounded as \(\hat{F}_S((S)) - \hat{F}_{S}(_{S}^*) \le G {(S)-_{S}^*} \le GD\) , where the first inequality uses \(G\) -Lipschitzness of \(\hat{F}_S\)  and the second the fact the both \((S)\)  and \(_{S}^*\)  lie in a ball of diameter \(D\) . Hence, attaining an excess empirical risk of \(GD\)  is trivial, and we now focus on deriving the other terms in the bounds.
Firstly, as discussed in [18]}, we consider \(G=D=1\) , since a simple reduction gives a factor of \(GD\)  for general \(G\)  and \(D\) . Furthermore, similar to [18]}, we show that the problem of TV-stable convex ERM is at least as hard as that of TV stable mean computation of a dataset with bounded mean - we state this reduction in prop:lowerboundreduction.
We now focus on showing accuracy lower bounds for \(\rho \) -TV-stable mean computation of dataset \(S\)  of size \(n\) , with mean \(\frac{M}{2}\le {\mu (S)} \le 2M\) . The accuracy, denoted by \(\alpha \) , is defined as \(\alpha ^2 = \mathbb {E}{(S)-\mu (S)}^2\) ,
\(\)  is a \(\rho \) -TV stable algorithm, and the expectation is taken over the algorithm's randomness.
The first part of thm:main-lower-bound follows thm:lowerbound11 which is based on a simple reduction argument.
This gives us that \(\alpha \ge \frac{1}{\rho n}\)  with \(M = \frac{1}{\rho n}\) . Plugging it in prop:lowerboundreduction, this gives us that excess empirical risk is lower bounded by \(\Omega {\frac{1}{\rho n}}\) .
Similarly, the second part follows from thm:lowerbound2 which gives us \(\alpha \ge \frac{1}{\sqrt{\rho n}}\)  with \(M=\frac{1}{\sqrt{\rho n}}\)  - the condition \(\alpha \le \frac{1}{4}\)  in the statement of thm:lowerbound2 can be absorbed in the trivial lower bound \(GD\) .
\(\Box \)

Discussion
In this work, we presented the TV stability framework for machine unlearning and instantiated it to develop unlearning algorithms for convex risk minimization problems. Currently, our results indicate two gaps, and motivate the following future directions.

Optimal TV-stable algorithm: Our upper and lower bound on excess empirical risk of TV stable algorithms don't match.
Hence, we either need to establish stronger lower bounds (arguably, more likely) or search for better algorithms.

Maximal coupling for unlearning: Our coupling procedure for unlearning for noisy-m-A-SGD is sub-optimal, in measure of its diagonal, by a \(\sqrt{T}\)  factor. A natural question is whether we can design an efficient maximal coupling.
We note that if efficiency were not a criteria, then this can be done - briefly, do a one step rejection sampling by computing the ratio of joint distribution iterates, if it fails, keep retraining, until the iterates generated is accepted by a rejection sampling.
However, in this case, the expected number of retrains can be shown to be one, and so is trivial. The challenge in this case is to give an efficient procedure when the first rejection sampling fails.

Beyond smooth convex functions:
The focus of this work was on smooth convex (loss) functions, but our techniques, and results for unlearning, extend to general non-convex functions.
However, a careful investigation of trade-offs between accuracy and unlearning efficiency, in classes of, say strongly-convex, non-smooth or even some non-convex functions, is an interesting future direction.


Acknowledgements
This research was supported in part by NSF BIGDATA award IIS-1838139 and NSF CAREER award IIS-1943251.

Additional related work
We survey the works on machine unlearning -
[27]} were one of the first papers to study the topic of machine unlearning. Their approach implements statistical query (SQ) algorithms by estimating the statistical queries using training data. Since the estimates are usually the mean of query evaluations computed on training data, unlearning is cheap, as we only need to subtract the evaluation on the deleted point.
[28]} studies this problem, with the goal to design systems to efficiently handle deletion requests.
Their approach, called SISA, is essentially a divide-and-conquer strategy, wherein the data is divided into disjoint sets, called shards, and a model on each shard is trained separately and aggregated. Furthermore, they do several check-pointing of states for each shard. In the average case, this provides a speedup of \(\frac{(R+1)S}{2}\)  for \(S\)  shards and \(R\)  checkpoints per shard, over retraining. They however give no guarantees on accuracy with this divide-and-conquer training method.
[2]} is another work which uses \((\epsilon , \delta )\) -differential privacy like guarantee. They study unlearning in generalized linear models, and propose a Newtons-step based method, leveraging connections with influence functions. Their computational cost is \(O(d^3)\)  computations for one unlearning. They, however give no guarantees on excess empirical risk achieved by the training method.
Finally, the work of [30]} studies batch unlearning in linear regression, with the goal to improve the computational cost of batch \(k\)  unlearning requests. Their method achieves a runtime of \(O(k^2d)\)  as opposed to \(O(kd^2)\)  for a naive approach. However, their notion of unlearning is again approximate, in the sense that model returned after unlearning is closest to the exact unlearning model among models in the \(d\)  dimensional subspace spanned by the to-be-deleted \(k\)  points. So it is easy to see that with larger \(k\) , the notion of approximation improves, which explains the \(k^2\)  term in the runtime as opposed to \(k\) .
Comparison with [1]}.
Our algorithm guarantees provable exact unlearning with probabilistic
runtime guarantees, whereas [1]} give algorithms with deterministic runtime and provide only an approximate \((\epsilon , \delta )\) -DP based unlearning guarantee – the \(\delta \)  can be interpreted as probability of the failure event in Monte-Carlo guarantees.
To handle these discrepancies when comparing, our stated runtime is the in-expecatation runtime. For a fixed runtime, we will look at regimes of \(\epsilon \)  and \(\delta \) , when the accuracy guarantee of [1]} is smaller than ours. We remind that a large \(\epsilon ,\delta \)  means a weaker unlearning criterion.
We will see that with same runtime, the accuracy of [1]} is smaller than ours in the regime when their unlearning parameters and hence the notion, is rather weak.
Considering the Lipschitz, smoothness parameters and diameter as constants, for smooth convex functions and \(k\)  edit requests, [1]} (Theorem 3.4) achieve an excess empirical risk of \(O{\frac{\sqrt{d}\sqrt{\log {1/\delta }}}{\epsilon n k}}^{2/5}\)  with an unlearning runtime of \(k^2\)  full-gradient computations. On the other hand, our algorithms achieve an an excess empirical risk of \(\min {\frac{1}{\sqrt{\rho n}}, {\frac{\sqrt{d}}{\rho n}}^{4/5}}\)  with \(\rho k\)  expected re-computations. Each re-computation takes \(m\cdot T\)  gradient computations where \(m\)  is the mini-batch size and \(T\)  the number of iterations. Therefore, in order to have the same runtime, we need \(\rho k m T= k^2n \iff \rho = \frac{kn}{mT}\) .
Firstly, note that as as long as \(d \le (\rho n)^{3/4}\) , noisy-m-A-SGD has smaller excess empirical risk than sub-sample-GD - this are the two regimes of interest.
We now set \(m\)  and \(T\)  for both the algorithms:
for Algorithm REF , \(m=\frac{\rho n}{T}\)  and \(T= \sqrt{\rho n}\) . This gives us \(\rho = \frac{kn}{\rho n} \iff \rho = \sqrt{k}\) , however \(\rho \)  is the Total Variation distance and is at most 1. Hence in regime \(d\ge (\rho n)^{3/4}\) , our runtime is always smaller than [1]}: \(kn\)  as opposed to \(k^2n\)  gradient computations. Even with \(\rho =1\) , our excess empirical risk is \(\frac{1}{\sqrt{n}}\)  and the excess empirical risk of [1]} is smaller than ours when \({\frac{\sqrt{d}\sqrt{\log {1/\delta }}}{\epsilon n k}}^{2/5} \lesssim \frac{1}{\sqrt{n}} \iff \frac{\epsilon }{\sqrt{\log {1/\delta }}} \gtrsim \frac{\sqrt{d}n^{1/4}}{k}\) .
In the second regime \(d< (\rho n)^{3/4}\) , we use Algorithm REF , wherein we have \(mT = \frac{(\rho n)^2}{d}\) . This gives us \(\rho = \frac{kn d}{(\rho n)^2} \iff \rho = \sqrt{\frac{kd}{n}}\) , and our excess empirical risk is \(O{{\frac{\sqrt{d}}{\rho n}}^{4/5}} = O{\frac{1}{(nk)^{2/5}}}\) .
Therefore, excess empirical risk of [1]} is smaller than ours when \({\frac{\sqrt{d}\sqrt{\log {1/\delta }}}{\epsilon n k}}^{2/5} \lesssim \frac{1}{(nk)^{2/5}} \iff \frac{\epsilon }{\sqrt{\log {1/\delta }}} \gtrsim d\) .
We therefore have that unless \(k\)  is very large, the accuracy of [1]} is smaller than ours when \(\epsilon \)  and \(\delta \) , take prohibitively large values which correspond to a weak notion of approximate unlearning.
We can similarly compare against Theorem 3.5 in [1]}, which will yield qualitatively similar conclusions.

Proofs for Section  REF
[Proof of Proposition REF ]
We first show that Algorithm REF  is \(\min {1,\rho }\) -TV stable for the aforementioned choice of number of iterations \(T\)  and mini-batch size \(m\) .
Consider neighbouring dataset \(S\)  and \(S^{\prime }\)  of \(n\)  points which differs in one sample, WLOG say the \(n^\text{th}\)  sample.
Let \((S):=\hat{}_S\)  and \((S^{\prime }):=\hat{}_{S^{\prime }}\)  denote the outputs of Algorithm REF  on \(S\)  and \(S^{\prime }\)  respectively.
Since in Algorithm REF , the randomness is only on indices, rather than actual data points, say that
\(S={1,2,\ldots ,n}\) .
Now we consider neighbouring dataset \(S^{\prime }\) , which contains \(n+1\)  or \(n-1\)  samples.
We will now consider the case when \(S^{\prime }\)  contains \(n-1\)  elements and the case with \(n+1\)  elements will follow analogously.
Let \(n\)  be the index present in \(S\)  but absent in \(S^{\prime }\) 
i.e. \(S^{\prime } = {1,2,\ldots , n-1}\) .
Let the sigma-algebra on these sets be the power sets of \(S\)  and \(S^{\prime }\)  respectively, denoted by \(\text{Pow}(S)\)  and \(\text{Pow}(S^{\prime })\)  respectively. Moreover, let \(\mu _{n,m}\)  denote the sub-sampling probability measure on \(n\)  points in \(S\)  i.e it sub-samples \(m\)  out of \(n\)  elements in \(S\)  uniformly randomly. Let \(\mu _{n,m}^{\otimes T}\)  denote the product measure of \(T\)  of \(\mu _{n,m}\) 's.
We similarly define \(\mu _{n-1,m}\)  and \(\mu _{n-1,m}^{\otimes T}\)  for \(S^{\prime }\) .
We first extend the sigma-algebra for the probability spaces so that the random variables \(\mu _{n,m}\)  and \(\mu _{n-1,m}\) , are defined on a common probability space.
For this, we will just add an event where the index \(n\)  can be sampled under \(\mu _{n-1,m}\)  with probability 0. We define \(\mu _{n,m}^{\prime }\)  as follows: for any set \(b \in \text{Pow}(S)\) ,
\(\mu _{n,m}^{\prime }(b) = {\left\lbrace \begin{array}{ll}\mu _{n-1,m}(b) & \text{ if } n \notin b\\ 0 & \text{otherwise} \end{array}\right.}\) .
We similarly extend the sigma algebra for the product space with measure
\(\mu _{n-1,m}^{\otimes T}\)  to get \(\mu _{n,m}^{^{\prime }\otimes T}\) .
Observe that for fixed initialization \(_0\)  and other parameters, Algorithm \((S)\)  and \((S^{\prime })\)  is the same (deterministic) map from \({\mathbf {b}} = (b_1,b_2,\ldots , b_T)\)  where \(b_j \in [n]^m\)  to \(\) .
They only differ because of different measures on the input space.
Hence total variation distance between \((S)\)  and \((S^{\prime })\)  is just the total variation distance between the push-forward measures \((S)_\# \mu _{n,m}^{\otimes T}\)  and \((S^{\prime })_\# \mu _{n,m}^{^{\prime }\otimes T}\)  which by using the fact that \((S)\equiv (S^{\prime }) \)  and data-processing inequality, is at most the total variation distance between \(\mu _{n,m}^{\otimes T}\)  and \(\mu _{n,m}^{^{\prime }\otimes T}\) .
Now the total variation distance can be bounded as,
\(\text{TV}((S), (S^{\prime })) & \le \text{TV}(\mu _{n,m}^{\otimes T}, \mu _{n,m}^{^{\prime }\otimes T}) =\sup _{{\mathbf {b}} \in \text{Pow}([n]^m)^T)} {\mu _{n,m}^{\otimes T}({\mathbf {b}}) - \mu _{n,m}^{^{\prime }\otimes T}({\mathbf {b}})}\\& = \mu _{n,m}^{\otimes T}\left({\mathbf {b}} \text{ such that at least one }b_j \text{ contains }n \right)\\& \le T \mu _{n,m}(b_1 \text{ contains } n) = \frac{Tm}{n}\) 

where the inequality follows using a union bound.
A similar argument works when \(S^{\prime }\)  is an neighbouring dataset of \(n+1\)  elements, yielding a total variation bound of \(\frac{Tm}{n+1} \le \frac{Tm}{n}\) .
Taking a uniform bound over all neighbouring datasets \(S^{\prime }\) , we get that \(\sup _{\Delta (S,S^{\prime })=1}\text{TV}((S),(S^{\prime })) \le \frac{Tm}{n}\) .
By definition of TV distance, we trivially have that \(\sup _{\Delta (S,S^{\prime })=1}\text{TV}((S),(S^{\prime }))\le 1\) . Therefore, setting \(m = \frac{\rho n}{T}\) , we get the desired result that the output of Algorithm REF  is \(\min {\rho ,1}\) -TV stable.
We now proceed to the accuracy guarantee which follows directly by analysis of SGD. We first show that the sub-sampling procedure produces unbiased gradients and bound its variance. For a fixed model \(\) , we have that
\(\mathbb {E}_b\frac{\sum _{j \in b} \nabla f(, z_j)}{m} = \sum _{\binom{n}{m} \text{ choices for }b}\frac{ \sum _{j \in b}\nabla f(, z_j)}{m \binom{n}{m} } = \frac{\binom{n-1}{m-1}}{m \binom{n}{m}} \sum _{j=1}^n \nabla f(, z_j) = \frac{\sum _{j=1}^n \nabla f(, z_j)}{n}\) 
where in the second equality, we use the observation that every \(_j\)  appears in exactly \(\binom{n-1}{m-1}\)  terms over all choices for \(b\) . We now bound its variance, denoted by a \(^2\)  by direct computation.
\(^2 &= \mathbb {E}_b {\frac{\sum _{j \in b} \nabla f(, z_j)}{m} - \mathbb {E}_b \left[\frac{\sum _{j \in b} \nabla f(, z_j)}{m} \right]}^2 \\&= \mathbb {E}_b {\frac{\sum _{j \in b} \nabla f(, z_j)}{m}}^2 - { \mathbb {E}_b \left[\frac{\sum _{j \in b} \nabla f(, z_j)}{m} \right]}^2 \\& = \sum _{{n \atopwithdelims ()m} \text{ choices for }b} \frac{1}{{n \atopwithdelims ()m}} \frac{1}{m^2} {\sum _{j \in b} \nabla f(, z_j)}^2 - {\frac{\sum _{j=1}^n \nabla f(, z_j)}{n}}^2\) 
In the first term, expanding the square and summing over all choices of \(b\) , we get exactly \({n-1 \atopwithdelims ()m-1}\)  terms of the form \({\nabla f(,_j)}^2\)  for \(j=1\)  to \(n\) , and \({n-2 \atopwithdelims ()m-2}\)  cross terms of the form \({\nabla f(,_i)}{\nabla f(,_j)}\)  for \(i \ne j\) , \(i, j=1\)  to \(n\) . Similarly, expanding the second term produces both these kind of terms. Accumulating the coefficients of all the terms, we get
\(^2 &=\mathbb {E}_b {\frac{\sum _{j \in b} \nabla f(, z_j)}{m} - \mathbb {E}_b \left[\frac{\sum _{j \in b} \nabla f(, z_j)}{m} \right]}^2\\&= {\frac{{n-1 \atopwithdelims ()m-1}}{m^2 {n \atopwithdelims ()m}} - \frac{1}{n^2}} \sum _{j=1}^n {\nabla f(,_j)}^2 + {\frac{{n-2 \atopwithdelims ()m-2}}{m^2 {n \atopwithdelims ()m}} - \frac{1}{n^2}}\sum _{i,j=1, i\ne j}^n{\nabla f(,_i)}{\nabla f(,_j)} \\& \le {\frac{1}{mn} -\frac{1}{n^2}}nG^2 + {\frac{m-1}{nm(n-1)} -\frac{1}{n^2}}\sum _{i,j=1, i\ne j}^n{\nabla f(,_i)}{\nabla f(,_j)} \\& \le {\frac{1}{m} -\frac{1}{n}}G^2 + {\frac{m-1}{nm(n-1)} -\frac{1}{n^2}}n(n-1)G^2\\& = {\frac{1}{m} -\frac{1}{n}}G^2 + {\frac{m-1}{m} -\frac{(n-1)}{n}}G^2 \\& = {\frac{1}{m} -\frac{1}{n}}G^2 + {\frac{1}{n}-\frac{1}{m}}G^2\\& = 2{\frac{1}{m} -\frac{1}{n}}G^2 \le \frac{2G^2}{m}\) 

where in the first inequality we used Cauchy-Schwartz inequality, and the fact the \(G\) -Lipschitzness implies the gradient norms are bounded by \(G\) . Finally, in the second last equality and the last inequality we used the fact that \(m\le n\) .
Since the sub-sampled gradients are unbiased, we can use the convergence guarantee of SGD on smooth convex function (see Theorem 4.1 in [41]}) which when using step size \(\eta \le \frac{1}{L}\)  gives us
\({ \hat{F}_S(\hat{}_S) - \hat{F}_S(^*_S)} \le O{\frac{ \eta ^2}{(1-\eta L)} + \frac{D^2}{\eta T}}\) 
Using step size \(\eta \le \frac{1}{2L}\) , the right hand side simplifies to \(2\eta ^2 + \frac{D^2}{\eta T} \le \frac{4G^2\eta }{m} + \frac{D^2}{\eta T} = \frac{4G^2T\eta }{\rho n} + \frac{D^2}{\eta T}\) , where in the last equality, we substituted \(m = \frac{\rho n}{T}\)  to ensure \(\rho \)  TV-stability. Balancing the trade off in \(\eta \)  gives us \(\eta = \frac{D\sqrt{\rho n}}{GT}\) . Therefore setting \(\eta = \min {\frac{1}{2L}, \frac{D\sqrt{\rho n}}{GT}}\)  gives us
\({ \hat{F}_S(\hat{}_S) - \hat{F}_S(^*_S)} \le O{\frac{GD}{\sqrt{\rho n}} + \frac{D^2L}{T}}\) 

Setting \(T = \frac{DL\sqrt{\rho n}}{G}\)  achieves the claimed result.
[Proof of Proposition REF ]
We first prove the stability guarantee. For this, we use the Rènyi-divergence based analysis used in differential privacy literature.
Let \(P\)  and \(Q\)  be probability distributions such that \(P\)  is absolutely continuous with respect to \(Q\)  and have densities \(\phi _P\)  and \(\phi _Q\) , respectively.
For \(\alpha \in (1,\infty )\) , the \(\alpha \)  Rènyi-divergence between \(P\)  and \(Q\)  is defined as follows [42]}:
\(D_\alpha (P\Vert Q) = \frac{1}{\alpha -1} \ln { \int f_P(x)^\alpha f_Q(x)^{1-\alpha } dx}\) 
Consider two neighbouring datasets \(S={_j}_{j}\)  and \(S^{\prime }={_j^{\prime }}_j\)  such that \(\Delta (S,S^{\prime })=1\) , and let \({b_t^{\prime }}_{t=1}^T\)  and \({_t^{\prime }}_{t=1}^T\)  denote the mini-batch indices and iterates of Algorithm REF  on dataset \(S^{\prime }\)  respectively. We look at iteration \(t\) , and fix all the randomness before \(t\)  i.e. fix \(_t \)  (and \(_{t}^{\prime }\) ), as well as randomness in sub-sampling mini-batch indices i.e. fix \(b_t\) . The \(\alpha \) -Rènyi Divergence between \(_{t+1}\)  and \(_{t+1}^{\prime }\)  can be bounded as,
\(D_\alpha (_{t+1} \Vert _{t+1}^{\prime }) &= D_\alpha {{ \mathring{}_t - \eta { \frac{\sum _{j \in b_t} \nabla f( \mathring{}_t,_j)}{m} +\theta _t}} \Big \Vert { \mathring{}_t - \eta { \frac{\sum _{j \in b_t} \nabla f( \mathring{}_t,_j^{\prime })}{m} +\theta _t}}}\\& \le D_\alpha { \mathring{}_t - \eta { \frac{\sum _{j \in b_t} \nabla f( \mathring{}_t,_j)}{m} +\theta _t} \Big \Vert  \mathring{}_t - \eta { \frac{\sum _{j \in b_t} \nabla f( \mathring{}_t,_j^{\prime })}{m} +\theta _t}}\\& \le D_\alpha { \frac{\sum _{j \in b_t} \nabla f( \mathring{}_t,_j)}{m} +\theta _t \Big \Vert \frac{\sum _{j \in b_t} \nabla f( \mathring{}_t,_j^{\prime })}{m} +\theta _t}\le \frac{2\alpha G^2}{m^2 \sigma ^2}\) 


where in the first and second inequality, we used post-processing property of Rènyi divergence, and in the last inequality, we use the fact that datasets \(S\)  and \(S^{\prime }\)  differ in at most one sample, therefore \({\frac{\sum _{j \in b_t} \nabla f( \mathring{}_t,_j)}{m} - \frac{\sum _{j \in b_t} \nabla f( \mathring{}_t,_j^{\prime })}{m} }^2 \le \frac{4G^2}{m^2}\) . Hence the divergence is between two multivariate Gaussians of same variance and with the square of the separation of their means at most \(\frac{2\alpha G^2}{m^2\sigma ^2}\) . Therefore, the inequality follows by using the formula for Rènyi divergence between two such multivariate Gaussians.
We now unfix \(b_t\) , and use the fact the \(b_t\)  is a uniform sample of \(m\)  out of \(n\)  (or \(n-1\)  or \(n+1\) ) indices. By privacy amplification by sub-sampling result in [43]}, for \(\alpha \le 2\) , we will argue that the Rènyi divergence upper bound amplifies to
\(\frac{32\alpha G^2}{n^2 \sigma ^2}\) . There are certain subtleties about the application of this result, so we explain, as follows. The first is that Theorem 9 stated in [43]},
when considering \(\alpha \le 2\) , the right hand side
simplifies as \( \frac{1}{\alpha -1}\log {1+\frac{m^2}{n^2}\frac{\alpha (\alpha -1)}{2}4{\exp {\frac{8G^2}{m^2\sigma ^2}}-1}} \le \frac{2\alpha m^2}{n^2}{\exp {\frac{8G^2}{m^2\sigma ^2}}-1} \le \frac{32\alpha G^2}{n^2 \sigma ^2}\)  where the last inequality use the numeric inequality \(\exp {x} \le 1+2x\)  when \(x\le 1.256\) ; this means that we need the following condition \(\frac{8G^2}{m^2\sigma ^2} \le 1.256\)  - we will revisit this condition later.
The second point is that Theorem 9 in [43]}
holds integer \(\alpha \ge 2\) , which only leaves us with \(\alpha =2\) .
In the subsequent part of the proof, we will need to take \(\alpha \rightarrow 1\) . This discrepancy can be accounted for by using the fact the \(\alpha \) -Rènyi Divergence is non-decreasing for \(\alpha \in [0,\infty ]\)  (see Theorem 3 in [46]}). Therefore the result holds for all \(\alpha \le 2\) , and we can replace the upper bound to be \(\frac{64G^2}{n^2\sigma ^2}\) 
The third and final point is that even though the amplification result in [43]} is established under the neighbouring relation that one point is replaced between datasets, it can be shown that the same result holds (perhaps upto constants) when the neighbouring relation is add/delete one data-point; see Lemma 3, [48]} for example.
We now use adaptive sequential composition property of Rènyi divergence (Proposition 1 in [49]}) which linearly accumulates the divergence across iterations, yielding that the Rènyi divergence between the iterates \((_1,_2,\ldots ,_T)\)  and \((_1^{\prime },_2^{\prime },\ldots ,_T^{\prime })\)  is bounded as, \(D_\alpha ((_1,_2,\ldots ,_T) \Vert (_1^{\prime },_2^{\prime },\ldots ,_T^{\prime })) \le \frac{64T G^2}{n^2 \sigma ^2}\) . An application of data-processing inequality gives us the same upper bound on the Rènyi divergence between the final iterates \(\hat{}_S\)  and \(\hat{}_S^{\prime }\) .
Moreover, this holds uniformly over all neighbouring datasets \(S^{\prime }\) . We now use the result that \(\lim _{\alpha \rightarrow 1}D_\alpha (\hat{}_S \Vert \hat{}_S^{\prime }) = D_{\text{KL}}(\hat{}_S \Vert \hat{}_S^{\prime })\)  where \(D_{\text{KL}}\)  denotes the KL-divergence (see Theorem 5 in [46]}). Hence we get that \( D_{\text{KL}}(\hat{}_S \Vert \hat{}_S^{\prime }) \le \frac{64TG^2}{n^2 \sigma ^2}\) . Finally, we use Pinsker's inequality to further lower bound the left hand side by total variation distance, which yields \(\text{TV}(\hat{}_S \Vert \hat{}_S^{\prime }) \le \sqrt{\frac{ D_{\text{KL}}(\hat{}_S \Vert \hat{}_S^{\prime })}{2}} \le \frac{8\sqrt{T}G}{n\sigma }\) . As remarked before, this is a uniform bound over all neighbouring datasets. Finally, as before, we trivially have that \(\text{TV}(\hat{}_S \Vert \hat{}_S^{\prime }) \le 1\) ; therefore setting \(\sigma = \frac{8\sqrt{T}G}{n\rho }\)  gives us that the algorithm's output is \(\min {\rho ,1}\)  TV-stable.
We now proceed to the accuracy guarantee. This follows simply by guarantee of Accelerated SGD on smooth convex functions. We have already shown in prop:upperboundsubsample that the gradients computed by sub-sampling are unbiased and its variance bounded by \(\frac{2G^2}{m}\) . The mean-zero Gaussian noise added preserves unbiasedness but the variance is bounded as,
\(^2 = {{\frac{\sum _{j \in b_t} \nabla f( \mathring{}_t,_j)}{m} +\theta _t -\nabla \hat{F}_S( \mathring{}_t) }}^2 &={{\frac{\sum _{j \in b_t} \nabla f( \mathring{}_t,_j)}{m} -\nabla \hat{F}_S( \mathring{}_t) }}^2 + {{\theta _t}}^2 \\& \le \frac{2G^2}{m} + \sigma ^2d\) 
We now use Theorem 2 from [51]} - they use notation \({\beta _t}_t\)  and \({\gamma _t}_t\)  for the step size schedule of Accelerated SGD and set \(\beta _t = \frac{t+1}{2}\)  and \(\gamma _t = \frac{t+1}{2}\gamma \) . Even though the updates of their A-SGD seem different than us, it can be verified that they are the same with \(\alpha _t = \beta _{t+1}(1-\beta _{t}^{-1}) = \frac{1-t}{t+2}\)  with \(\alpha _0=0\)  and \(\eta = \gamma \) .
Finally, using step-size \(\eta \le \frac{1}{2L}\) , and appealing to Theorem 2 in [51]}, we get,
\({ \hat{F}(\hat{}_S) - \hat{F}(^*)} \le O{ T\eta ^2 + \frac{D^2}{\eta T^2}} = O{\eta T {\frac{2G^2}{m} +\sigma ^2d} + \frac{D^2}{\eta T^2}}\) 
Let \(\tilde{G}^2 = \frac{2G^2}{m} +\sigma ^2d\) , balancing the trade-off in \(\eta \)  gives us \(\eta = \frac{D}{\tilde{G} T^{3/2}}\) . Therefore, setting \(\eta = \min {\frac{1}{2L},\frac{D}{\tilde{G}T^{3/2}}}\)  gives us
\({ \hat{F}(\hat{}_S) - \hat{F}(^*)}&\le O{\frac{LD^2}{T^2} + \frac{\tilde{G}D}{\sqrt{T}}} \le O{\frac{LD^2}{T^2} + \frac{GD}{\sqrt{Tm}} + \frac{\sigma \sqrt{d} D}{\sqrt{T}}} \\&\le O{\frac{LD^2}{T^2} + \frac{GD}{\sqrt{Tm}} + \frac{G D\sqrt{d}}{n \rho }}\) 
Finally, note that when using the amplification lemma, we arrived at the condition \(\frac{8G^2}{m^2\sigma ^2} \le 1.256\) . Substituting \(\sigma = \frac{8\sqrt{T}G}{n\rho }\) , this reduces to \(\frac{(n\rho )^2}{8m^2T} \le 1.256 \iff T \ge \frac{(n\rho )^2}{16m^2}\) .
[Proof of corr: upper-bound-m-A-sgd]
We start with the result in Proposition REF , and balance the two trade-offs:
the first between the terms \(\frac{GD}{\sqrt{mT}}\)  and \(\frac{GD\sqrt{d}}{\rho n}\) , and the second between \(\frac{GD\sqrt{d}}{\rho n}\)  and \(\frac{LD^2}{T^2}\) .
Note that as long as \(\frac{GD}{\sqrt{mT}} \ge \frac{LD^2}{T^2} \iff m \le \frac{T^3 G^2}{(LD)^2}\) , the second term is larger than the first.
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}\) .
Similarly, optimizing the trade-off between the first and third term gives us
\(\frac{GD\sqrt{d}}{\rho n} = \frac{LD^2}{T^2} \iff T = \sqrt{\frac{LD(\rho n)}{G\sqrt{d}}}\) . Hence setting \(T = \max {\frac{(\rho n)^2}{md},\sqrt{\frac{LD(\rho n)}{G\sqrt{d}}} }\)  yields an expected excess empirical risk of \(O{\frac{GD\sqrt{d}}{n \rho }}\) .
We now look at the given condition \(T\ge \frac{(n\rho )^2}{16m^2}\)  given in prop:upper-bound-noisy-m-a-sgd.
We have set \(T = \max {\frac{(\rho n)^2}{md},\sqrt{\frac{LD(\rho n)}{G\sqrt{d}}} }\) , there we need to ensure that \(\frac{(\rho n)^2}{md} \ge \frac{(\rho n)^2}{16m^2} \iff m \ge \frac{d}{16}\) , as well as \(\sqrt{\frac{LD(\rho n)}{G\sqrt{d}}} \ge \frac{(\rho n)^2}{16m^2} \iff m \ge \frac{1}{4}{\frac{(\rho n)^3 G\sqrt{d}}{LD}}^{1/4}\)  - this recovers the condition \(m \ge \min {\frac{d}{16}, \frac{1}{4}{\frac{(\rho n)^3 G\sqrt{d}}{LD}}^{1/4}}\)  in the Proposition statement.
Combining all the above arguments, we get that for any \(m \ge \min {\frac{d}{16}, \frac{1}{4}{\frac{(\rho n)^3 G\sqrt{d}}{LD}}^{1/4}}\)  , setting \(T = \max {\frac{(\rho n)^2}{md},\sqrt{\frac{LD(\rho n)}{G\sqrt{d}}}}\) , yields an expected excess empirical risk of \(O{\frac{GD\sqrt{d}}{n \rho }}\) .
Note that in the above proof, if we use the stronger variance bound of \(2L^2{\frac{1}{m}-\frac{1}{n}}\)  from sub-sampling (derived in the proof of Proposition REF ), we get that when doing full-gradient descent, the variance, as expected is zero, which yields a running time of \(T = \sqrt{\frac{LD\rho n}{G\sqrt{d}}}\) .
There exists neighbouring datasets \(S\)  and \(S^{\prime }\)  of \(n\)  points, and smooth \(G\) -Lipshcitz convex functions \(f\)  and constraint set \(\)  such that the total variation distance between iterates produced by Algorithm REF  run on datasets \(S\)  and \(S^{\prime }\) , denoted by \({_1,_2,\ldots , _T}\)  and \({_1^{\prime },_2^{\prime },\ldots , _T^{\prime }}\)  respectively, is bounded as \(\text{TV}({_1,_2,\ldots , _T}, {_1^{\prime },_2^{\prime },\ldots , _T^{\prime }}) \ge \min {\Omega {\frac{G\sqrt{T}}{n\sigma }},1}\) .
[Proof of prop:tv-stability-noisy-sgd]
We first prove this without projection - let the constraint set \(=^d\) , and so the projection \(\)  is the identity map. Also, for simplicity, let the initial model be 0. Consider data sets \(S\)  and \(S^{\prime }\)  such that all points are 0 but the \(n^{\text{th}}\)  differing point. Let the \(n^{\text{th}}\)  point of \(S\)  be \(-G_1\)  and that of \(S^{\prime }\)  be \(G_1\) , where \(_1\)  is the first canonical basis vector. Let the function \(f(,) = {}{}\) . The gradients are just data points \(\) , therefore gradients are 0 on all but the differing points, wherein in the differing point in dataset \(S\) , the gradient is a constant \(-G_1\)  and for dataset \(S^{\prime }\) , it is \(G_1\) .
Consider the map \(\Psi : (x_1,x_2,\ldots , x_T) \rightarrow x_T\) ; using data processing inequality and this map, we have that
\(\text{TV}({_1,_2,\ldots , _T}, {_1^{\prime },_2^{\prime },\ldots , _T^{\prime }}) \ge \text{TV}(\Psi {_1,_2,\ldots , _T}, \Psi {_1^{\prime },_2^{\prime },\ldots , _T^{\prime }}) = \text{TV}(_T, _T^{\prime })\) 

We now focus on bounding the total variation distance between the last iterates. Furthermore, by data-processing inequality, we can get rid of the step size scaling, and therefore can consider the last iterates as just the sum of all gradients. By simple calculations, we get that \(_T\)  is a mixture of multivariate Gaussians, all with variance \(T\sigma ^2 \)  but with varying means: \(G_1, 2G_1, \ldots , TG_1\) , similarly for \(_T^{\prime }\) . We denote the mixtures probabilities by \(\pi _i\)  where the \(i^{th}\)  conditional distribution, denoted by \(_T^i\)  and \(_T^{\prime i}\)  respectively, has means \(iG_1\)  and \(-iG_1\)  respectively. Also, we denote the conditional probability densities of the \(i^{th}\)  distribution by \(\phi _S^i()\)  and \(\phi _{S^{\prime }}^i()\)  respectively.
We will show that the total variation between these mixtures is expected total variation distance between the mixture components. This follows due the symmetry between these two mixtures, which implies that the set that achieves the total variation distance is \({: _1 \ge 0}\) .
We can therefore write the total variation distance as,
\(\text{TV}(_T \Vert _T^{\prime }) &= \frac{1}{2} {\phi _S() - \phi _{S^{\prime }}()}_1 = \int _{_1\ge 0} \phi _S() -\phi _{S^{\prime }}()d= \int _{_1\ge 0} \sum _{i} \pi _i (\phi ^i_S() -\phi ^j_{S^{\prime }}())d\\& = \sum _{i} \pi _i \int _{\ge 0} (\phi ^i_S() -\phi ^j_{S^{\prime }}())d= \sum _{i}\pi _i \text{TV} (_T^i, _T^{\prime i}) \\& \gtrsim \sum _i \pi _i \frac{2Gi}{m \sqrt{T}\sigma } = \frac{2G}{m\sqrt{T}\sigma }{i} =\frac{2G\sqrt{T}}{n\sigma }\) 
where in the inequality, we use the fact that \(_T^i\)  and \(_T^{\prime i}\)  are Gaussians with means separated by \(2G\) , and variance being \(T\sigma ^2\)  and use the lower bound result on TV between high-dimensional Gaussians [53]}. Finally, in the last equality, we compute the Expected value of \(i\)  under the mixture distribution - recall that \(i\)  is a sum of \(T\)  Bernoulli random variables with bias \(\frac{m}{n}\) , the expectation of which is \(\frac{Tm}{n}\) .
We now argue why projection doesn't change the above claim. Note the with the projection, all the Gaussians in the mixture are truncated forming a discrete distributions at the boundary of the constraint set.
The probability mass on either sides of the (original) mean is unchanged. Hence \({: _1 \ge 0}\)  is still the witness set of total variation distance between the mixtures, and the total variation distance in both constrained/unconstrained cases is the same.
The same holds for the total variation between the corresponding mixture components.
These observations suffices for application of proof of the unconstrained case. Finally, since TV distance, by definition is upper bounded by 1 - this gives a trivial lower bound of 1, and hence the TV distance is lower bound by \(\min {\Omega {\frac{G\sqrt{T}}{n \rho }},1}\) .

Proofs for Section  REF
We introduce some notation and setup the roadmap.
In the start of the stream, we have a model trained on the initial dataset of \(n\)  samples. We then observe an insertion or deletion request. We enumerate the data points from 1 to \(n\) , and without loss of generality, assume that the \(n^\text{th}\)  sample is to be deleted, and the inserted sample has index \(n+1\) .
We want to show that the unlearning algorithm satisfies exact unlearning at every time point in the stream, and what suffices is to argue that this holds for one edit request, since by mathematical induction it then holds for the entire stream.
For one edit request, we will show the following:
1. unlearning (deletion/insertion) algorithm is a valid transport, and
2. the probability of recompute is small, and we will see that together these will imply, that it is a coupling, with large enough measure of the diagonal.
Let \(\mu _{n,m}\)  denote the sub-sampling probability measure
to sample \(m\)  out of \(n\)  elements uniformly randomly.
In the deletion and insertion algorithms, we replace some mini-batch indices in some iterations - let these operations be denoted by \(\textsc {Del}\)  and \(\textsc {Ins}\)  respectively.
To elaborate, \(\textsc {Del}\)  is a (deterministic) map from \({[n]^m}^T\)  to \({[n-1]^m}^T\)  and \(\textsc {Ins}\)  is a map from \({[n]^m}^T\)  to \({[n+1]^m}^T\) .
For an input \(b \in ([n]^m)^T\) , we have that \({\mathbf {b}} \sim \mu _{n,m}^{\otimes T}\) .
Furthermore, define \(\mu _{n,m}^{\text{del} \otimes T} := \textsc {Del}_\# \mu _{n,m}^{\otimes T}\)  and \(\mu _{n,m}^{\text{ins} \otimes T} := \textsc {Ins}_\# \mu _{n,m}^{\otimes T}\) .
An important observation is that in the unlearning Algorithm REF , the sub-sampled indices \({\mathbf {b}}\)  are drawn from a product distribution \( \mu _{n,m}^{\otimes T}\)  and in each iteration of alg:unlearn-subsample-gd or alg:unlearn-noisy-m-a-sgd, the maps \(\textsc {Del}\)  and \(\textsc {Ins}\)  act component-wise and symmetrically.
This implies that \(\textsc {Del}({\mathbf {b}}) = [\text{del}(b_1),\text{del}(b_2),\ldots , \text{del}(b_T)]\)  where \(\text{del}: [n]^m \rightarrow [n-1]^m\)  is the function which describes one iteration of the unlearning algorithm for handling mini-batch indices.
We similarly have function \(\text{ins}: [n]^m \rightarrow [n+1]^m\)  for insertion.
We finally define \(\mu _{n,m}^{\text{del}} := \text{del}_\# \mu _{n,m}\)  and \(\mu _{n,m}^{\text{ins}} := \text{ins}_\# \mu _{n,m}\) 
- these are the probability measures induced on the sub-sampling indices by deletion and insertion operations, respectively.
Unlearning for sub-sample-GD
We first show that \(\mu ^{\text{del}}_{n,m}\) , the probability distribution, induced at a given iteration during deletion, over mini-batch indices \(b \in [n]^m\)  is a transport.
[Deletion]

For any set \(b \in [n]^m\) , we have that \(\mu ^{\text{del}}_{n,m}(b) = \mu _{n-1,m}(b)\)
First note that if the verification is unsuccessful, then a recompute is triggered and therein at each iteration, we drawn \(b \sim \mu _{n-1,m}\) .
Therefore, \(\mu ^{\text{del}}_{n,m}(b)= \mu _{n-1,m}(b)\)  follows trivially.
We now argue for the other case.
The verification is successful if the deleted point was not present in any of iterations, i.e. at any iteration the sub-sample batch \(b_t\)  doesn't contain the deleted point \(\) .
The measure \(\mu _{n,m}^{\text{del}}\)  is therefore just the probability under the original sub-sampling measure \(\mu _{n,m}\)  conditioned on the event that \(\notin b\) .
We therefore have,
\(\mu _{n,m}^{\text{del}}(b) &= \mu _{n,m}(b|{\notin b}) =\frac{ \mu _{n,m}(b \cap {\notin b}) }{\mu _{n,m}({\notin b})}\) 

By direct computation, \(\mu _{n,m}({\notin b}) = 1-\mu _{n,m}({\in b}) = 1 - \frac{{n-1 \atopwithdelims ()m-1}}{{n \atopwithdelims ()m}} = 1- \frac{m}{n}\) . We now look at two choices for \(b\) . First suppose \(z \in b\) , then the numerator \(\mu _{n,m}(b \cap {\notin b}) = 0\) , which gives us that \( \mu _{n,m}^{\text{del}}(b) = 0 = \mu _{n-1,m}(b)\) . We now look at a \(b\)  such that \(z \notin b\) . We have,
\(\mu _{n,m}^{\text{del}}(b) &=\frac{ \mu _{n,m}(b) }{\mu _{n,m}({\notin b})}=\frac{1/{n \atopwithdelims ()m}}{1- m/n} = \frac{n}{n-m}\frac{(n-m)!m!}{n!} \\&=\frac{(n-m-1)!m!}{(n-1)!} = \frac{1}{{n-1 \atopwithdelims ()m}} = \mu _{n-1,m}(b)\) 
Similarly, for insertion, we show that
\(\mu ^{\text{ins}}_n\) , the probability distribution, induced at a given iteration during insertion, over mini-batch indices \(b \in [n]^m\) , is a is valid transport.
For any set \(b \in [n+1]^m\) , we have that \(\mu ^{\text{ins}}_{n,m}(b) = \mu _{n+1,m}(b)\) 

[Proof of claim:sub-sample-gd-insertion]
Let \(\nu \)  denote the uniform probability measure over \(n+1-m\)  elements.
Given \(b\) , we consider two cases based of whether last/inserted index \(n+1\)  lies in \(b\)  or not.
In the first case, we know that the outcome of \(\text{Bernoulli}(m/(n+1))\)  must have been 1 i.e. the iteration was selected. Furthermore, in that case, the inserted point would have replaced some other point not in \(b\)  - the total number of possibilities are \(n+1-m\) . Let \(E_i\)  be event that the inserted point replaced the \(i^{\text{th}}\)  data point, whose index we denote by \(s_i\) . Note that the events \(E_i^{\prime }s\)  are disjoint and the event \(b\)  is \(\cup _{i=1}^{n+1-m} E_i\) .
Furthermore, \(\mu ^{\text{ins}}_{n,m}(E_i) = \mu ^{\text{ins}}_{n,m}(\text{original subsample is } b\backslash {n+1}\cup {s_i} | {s_i} \text{ replaced}) \mu ^{\text{ins}}_{n,m}({s_i} \text{replaced})) = \mu _{n,m}(b\backslash {n+1}\cup {s_i} | {s_i}) \nu ({s_i}) = \frac{1}{{n \atopwithdelims ()m-1}}\frac{1}{(n+1-m)}\) .
We therefore have that
\(\mu ^{ins}_{n}(b) &=\frac{m}{n+1} \mu ^{ins}_{n}(\cup _{i=1}^{n+1-m} E_i) = \frac{m}{n+1}\sum _{i=1}^{n+1-m}\mu ^{ins}_{n}(E_i) = \frac{m}{n+1}\sum _{i=1}^{n+1-m}\frac{1}{{n \atopwithdelims ()m-1}}\frac{1}{(n+1-m)} \\& = \frac{m}{n+1}\frac{1}{{n \atopwithdelims ()m-1}} = \frac{m (m-1)! (n-(m-1))!}{(n+1) n!} = \frac{1}{{n+1 \atopwithdelims ()m}} = \mu _{n+1,m}(b)\) 

In the other case, we know that Bernoulli(\(m/(n+1)\) ) resulted in 0, so there is no replacement. Therefore, we have
\(\mu ^{ins}_{n}(b) = {1-\frac{m}{n+1}}\frac{1}{{n \atopwithdelims ()m}} = \frac{(n+1-m)(n-m)!m!}{n!(n+1)} = \frac{1}{{n+1 \atopwithdelims ()m}} = \mu _{n+1,m}(b)\) 
Coupling.
We formally describe the coupling constructed by the unlearning alg:unlearn-subsample-gd. We first the discuss deletion case - consider datasets \(S\)  and \(S^{\prime }\)  of sizes \(n\)  and \(n-1\)  respectively, and wlog assume that the last sample of \(S\)  differs.
We first sample \(\mathbf {b} = [b_1,b_2,\ldots , b_T] \sim \mu _{n,m}^{\otimes T}\) .
We set \(\mathbf {b}^{(1)}=\mathbf {b}\) .
For each \(j \in T\) ,
if \(n \in b_j\) , then sample \(b^{(2)}_j \sim \mu _{n-1,m}\) , otherwise set \(b^{(2)}_j = b_j\) .
This produces the coupled mini-batches \((\mathbf {b}^{(1)}, \mathbf {b}^{(2)})\)  for deletion.
For insertion, we have datasets \(S\)  and \(S^{\prime }\)  of sizes \(n\)  and \(n+1\)  respectively, and again assume that the last of point of \(S^{\prime }\)  differs. Sample \(\mathbf {b} = [b_1,b_2,\ldots , b_T] \sim \mu _{n,m}^{\otimes T}\) .
and set \(\mathbf {b}^{(1)}=\mathbf {b}\) .
Now sample \({c_{j}}_{j=1}^{T}\) , where \(c_{j} \sim \text{Bernoulli}{\frac{m}{n}}\) , if \(c_{j}=1\) , then sample uniformly a point in \(b^{(2)}_j\) , and replace it with \(n+1\) . Otherwise set \(b^{(2)}_j = b_j\) , which gives us the coupled mini-batches \((\mathbf {b}^{(1)}, \mathbf {b}^{(2)})\) .
It is easy to see that the above procedure is how Algorithm REF  handles insertions and deletions going from \(S\)  to \(S^{\prime }\) .
We first show that this is a valid coupling.
For the coupling described above, for any \(b\) ,

\( ¶{\mathbf {b}^{(1)} = b} = \mu _{n}^{\otimes T}(b)\)

\(¶{\mathbf {b}^{(2)} = b} = \mu _{n-1}^{\otimes T}(b) \ \ \text{(deletion)}\) ,   \(¶{\mathbf {b}^{(2)} = b} = \mu _{n+1}^{\otimes T}(b) \ \ \ \text{(insertion)}\)

[Proof of claim:subsample-gd-coupling]
Follows immediately from Claims REF  and REF  .
We now show that the probability of disagreement under the above coupling is upper bounded by \(k\)  times TV-stability parameter of alg:sub-sample-GD.
For the \(\rho \) -TV stable alg:sub-sample-GD,
under the coupling described above, the following holds
\(_{(\mathbf {b}^{(1)},\mathbf {b}^{(2)})}[\mathbf {b}^{(1)} \ne \mathbf {b}^{(2)}] \le \rho \) 
For deletion, we have,
\(\mathbb {P}_{(\mathbf {b}^{(1)}, \mathbf {b}^{(2)})}[\mathbf {b}^{(1)} \ne \mathbf {b}^{(2)}] = \mathbb {P}_{(\mathbf {b}^{(1)}, \mathbf {b}^{(2)})}[ \exists j \in [T]: \mathbf {b}^{(1)}_j \ne \mathbf {b}^{(2)}_j] = \mathbb {P}_{\mathbf {b}}[\exists j \in [T]:n \in \mathbf {b}_j]\le \frac{Tm}{n}\) 
For insertion, we have
\(\mathbb {P}_{(\mathbf {b}^{(1)}, \mathbf {b}^{(2)})}[\mathbf {b}^{(1)} \ne \mathbf {b}^{(2)}] &= _{{\mathbf {b}}, \mathbf {c}}[ \exists j \in [T] : c_j=1] \le \frac{Tm}{n}\) 
In Proposition REF , we showed that the total variation distance of the algorithm under change of one point is at most \(\frac{Tm}{n} =\rho \) , which completes the proof.
We are now ready to prove prop:unlearn-subsample-GD.
[Proof of prop:unlearn-subsample-GD]
The following argument is for deletion, but the insertion case follows similarly.
Consider dataset \(S\)  and \(S^{\prime }\)  of \(n\)  points and \(n-1\)  points respectively, differing in one sample.
As in the proof of prop:upperboundsubsample, we embed the randomness for Algorithm REF  executed on \(S\)  and \(S^{\prime }\)  into a common probability space.
Therefore, similar to the proof of prop:upperboundsubsample given the datasets (and other parameters),
Algorithm REF , \((S)\) 
is a deterministic map from sub-sampled indices \({\mathbf {b}} = (b_1,b_2,\ldots , T)\)  to the model: \((S) :{\mathbf {b}} \rightarrow \) , where \(b_j \in [n]^m\) , for both datasets.
Hence, what suffices is to show that the input probability measure \(\mu _{n,m}^{\otimes T}\)  is transported to the one that would have been produced on the current dataset \(S^{\prime }\)  i.e \(\mu _{n-1,m}^{\otimes T}\)  - this follows from claim:subsample-gd-coupling.
Hence it follows that the output generated by Algorithm REF  has the same measure as \((S)_\# \mu _{n-1,m}^{\otimes T}\) , which proves first part of the claim.
The probability of recompute, being at most \(\rho \) , for one edit, follows directly from REF . Finally, from remark:tvuptok, for \(k\)  edits, and the assumption the number of samples throughout the stream is between \(n/2\)  and \(2n\) , the recompute probability is at most \(2k\rho \) .

Unlearning for noisy-m-A-SGD
Coupling mini-batches
In this section, we show that Algorithm REF  transports sub-sampling probability measures while handling edit requests.
We remind that \(\mu ^{\text{del}}_{n,m}\)  denotes the probability measure induced on the sub-sampled indices by the deletion procedure, in any iteration.
We show that, for any mini-batch, the probability mass of the mini-batched indices under \(\mu ^{\text{del}}_{n,m}\)  is same as that under the sub-sampling measure \(\mu _{n-1,m}\) .
For any set \(b \in [n-1]^m\) , we have that \(\mu ^{\text{del}}_{n,m}(b) = \mu _{n-1,m}(b)\) 

[Proof of claim:sgddeletion]
Firstly, note that deletion uses additional randomness which is used to uniformly sample one element from \(n-(m-1)\)  elements - let \(\nu \)  denote the uniform probability measure on \(n-(m-1)\)  elements.
Let \(E\)  be the event that the \(n^{th}\)  was sub-sampled originally, and therefore replaced upon verification. By direct computation \(\mu _{n,m}(E) = \frac{m}{n}\) . We can therefore write \(\mu ^{\text{del}}_{n,m}(b)\)  as follows
\(\mu ^{\text{del}}_{n,m}(b) &= \mu ^{\text{del}}_{n,m}(b|E) \mu _{n,m}(E) + \mu ^{\text{del}}_{n,m}(b|E^c) \mu _{n,m}(E^c) \\\) 

Under event \(E\) , we have the deleted index was replaced. But it can be any element of \(b\)  that arised out of this replacement. Hence we decompose the event \(b|E\)  into events \(E_i\) 's, where \(E_i\)  corresponds to the event that \(b_i\)  was replaced. We have that \(b|E = \cup _{i=1}^m E_i\) , and furthermore, due to the uniform measure, \(\mu ^{\text{del}}_{n,m}(E_i) = \mu ^{\text{del}}_{n,m}(E_j) \forall i,j\) . Note that in the event \(E_i\) , we require that the original sub-sampling measure on \(n\)  points \(\mu _{n,m}\)  to have produced the set \(b\backslash b_i \cup {n}\)  and then a uniform \(b_i\)  is drawn upon replacement.
Therefore, \(\mu ^{\text{del}}_{n,m}(E_i) = \mu _{n,m}(b \backslash b_i \cup {n}) \nu (b_i) = \frac{1}{{n-1 \atopwithdelims ()m-1}} \frac{1}{n-1-(m-1)}\) . Similarly, when the event \(E^c\)  occurs, probability of outputting \(b\)  corresponds to the event when \(b\)  was generated using the original sub-sampling measure \(\mu _m\)  (and no additional randomness used upon verification). Therefore, we get \(\mu ^{\text{del}}_{n,m}(b|E^c) = \mu _{n,m}(b | E^c) = \frac{1}{{n-1 \atopwithdelims ()m}}\) . Plugging these in, and with simple calculations, we have
\(\mu ^{\text{del}}_{n,m}(b) &= \sum _{i=1}^m\mu ^{\text{del}}_{n,m}(E_i) \mu _{n,m}(E_i) + \mu ^{\text{del}}_{n,m}(b|E^c) \mu _{n,m}(E) \\&= \sum _{i=1}^m\mu _{n,m}(b \backslash b_i \cup {n}) \nu (b_i) \frac{m}{n} + \frac{1}{{n-1 \atopwithdelims ()m}} {1-\frac{m}{n}} \\& =\frac{m}{{n-1 \atopwithdelims ()m-1}} \frac{1}{n-1-(m-1)}\frac{m}{n}+ \frac{1}{{n-1 \atopwithdelims ()m}} {1-\frac{m}{n}}\\& = \frac{1}{{n-1 \atopwithdelims ()m}} + \frac{m}{n}{\frac{m}{{n-1 \atopwithdelims ()m-1}(n-1-(m-1))} - \frac{1}{{n-1 \atopwithdelims ()m}}} \\& = \frac{1}{{n-1 \atopwithdelims ()m}} + \frac{m}{n}{\frac{m(m-1)!(n-1-(m-1))!}{(n-1)!(n-1-(m-1))} - \frac{1}{{n-1 \atopwithdelims ()m}}} \\& = \frac{1}{{n-1 \atopwithdelims ()m}} + \frac{m}{n}{\frac{1}{{n-1 \atopwithdelims ()m}} - \frac{1}{{n-1 \atopwithdelims ()m}}} = \frac{1}{{n-1 \atopwithdelims ()m}} = \mu _{n-1,m}(b) \) 
Similarly, for insertion, we now show that the probability mass of any mini-batch under \(\mu ^{\text{ins}}_{n,m}\) , the probability measure induced by insertion on the \(n+1\)  data points, is same as that under \(\mu _{n+1,m}\) .
For any set \(b \in [n+1]^m\) , we have that \(\mu ^{\text{ins}}_{n,m}(b) = \mu _{n+1,m}(b)\) 

[Proof of claim:sgdinsetion]
Same as that of claim:sub-sample-gd-insertion.

Lemmas for reflection coupling
We state and prove some results about reflection mapping and couplings.


Let \(P\)  and \(Q\)  be probability distributions over \(^d\) .
Let \(\psi :^d\rightarrow ^d\)  be a bijection such that \(\phi _P(\psi (x))= \phi _Q(x)\) , \(\phi _P(\psi ^{-1}(x)) = \phi _Q(x)\)  and \({\det \left(\frac{d\psi (x)}{dx}\right)}=1\) , where \(\frac{d\psi (x)}{dx}\)  is the Jacobian of the multivariate map \(\psi \) .
Let \(x \sim P\)  be a sample from \(P\) . Let \(y=x\)  if \(\text{Unif}(0,1)\le \frac{\phi _Q(x)}{\phi _P(x)}\) , otherwise \(y=\psi (x)\) . Then \((x,y)\)  is a maximal coupling of \(P\)  and \(Q\) .


We first show that \(y\)  is a sample from \(Q\) . Let \(E\)  be an event in the range of \(Q\) . Let accept be the event when \(u \sim \text{Unif}(0,1)\) , \(u \le \frac{\phi _Q(x)}{\phi _P(x)}\) . We have,
\(¶{y \in E} & = ¶{y\in E, \text{accept}} + ¶{y\in E, \text{reject}} \\& = \mathbb {E}_{x,u}\left[\mathbb {1}{x \in E} \mathbb {1}{u \le \frac{\phi _Q(x)}{\phi _P(x)}}\right] + \mathbb {E}_{x,u}\left[\mathbb {1}{\psi (x) \in E} \mathbb {1}{u > \frac{\phi _Q(x)}{\phi _P(x)}}\right] \\& = \mathbb {E}_{x}\left[\mathbb {1}{x \in E} ¶{{u \le \frac{\phi _Q(x)}{\phi _P(x)}} \Big \vert x}\right] + \mathbb {E}_{x}\left[\mathbb {1}{\psi (x) \in E} ¶{{u > \frac{\phi _Q(x)}{\phi _P(x)}}\Big \vert x}\right] \\& = \int _{^d}\mathbb {1}{x \in E} \min {1, \frac{\phi _Q(x)}{\phi _P(x)}}\phi _P(x)dx \\&+ \int _{^d}\mathbb {1}{\psi (x) \in E} {1 - \min {1, \frac{\phi _Q(x)}{\phi _P(x)}}}\phi _P(x)dx \\& = \int _{^d}\mathbb {1}{x \in E} \min {\phi _P(x), \phi _Q(x)}dx + \int _{^d}\mathbb {1}{\psi (x) \in E} \max {0,\phi _P(x) - \phi _Q(x)}dx\) 

For the second term, we now do change of variable - let \(v = \phi (x)\)  - using the given properties of \(\psi \) , we have \(\phi _P(x) = \phi _P(\psi ^{-1}(v)) = \phi _Q(v)\)  and \(\phi _Q(x) = \phi _P(v)\) . Furthermore \(dv = {\det \left(\frac{d\psi (x)}{dx}\right)}dx = dx\) . Finally, we are integrating over \(^d\) , and since \(\phi \)  is a bijection, it can flip the limits of some of the coordinates, however, that is taken into account with using the absolute value of the determinant of the Jacobian.
The second term therefore becomes \(\int _{^d}\mathbb {1}{v \in E} \max {0,\phi _Q(v) - \phi _P(v)}dv\) . We now combine the integrands of both the terms, and substitute \(v=x\)  as the variable in the second term. This gives us,
\(¶{y \in E} & = \int _{^d} \mathbb {1}{x \in E} {\min {\phi _P(x), \phi _Q(x)} + \max {0,\phi _Q(x) - \phi _P(x)}}dx\) 

Note that for a fixed \(x\) , if \(\phi _P(x) \le \phi _Q(x)\) , the integrand becomes \( \mathbb {1}{x \in E} {\phi _P(x) + \phi _Q(x) - \phi _P(x)} = \mathbb {1}{x \in E} Q(x)\) . On the other hand, if \(\phi _P(x) > \phi _Q(x)\) , the integrand becomes \(\mathbb {1}{x \in E}\phi _Q(x)\) . Hence, for all cases, we get that,
\( ¶{y \in E} = \int _{^d} \mathbb {1}{x \in E} \phi _Q(x)dx = Q(E)\) 
We now show that it is a maximal coupling i.e. the probability of accept is \(1-\text{TV}(P,Q)\) . We have,
\(¶{\text{accept}} &= \mathbb {E}_{x,u}\left[\mathbb {1}{u \le \frac{\phi _Q(x)}{\phi _P(x)}}\right] = \int _{^d} \min {1,\frac{\phi _Q(x)}{\phi _P(x)}}\phi _P(x)dx \\& = \int _{^d} \min {\phi _P(x),\phi _Q(x)}dx = 1-\text{TV}(P,Q) \) 
Let \(P\)  and \(Q\)  be two isotropic probability distributions over \(^d\)  with means \(\mu _P\)  and \(\mu _Q\)  such that for any vectors \(,\) , \(\phi _P() = \phi _Q()\)  if \({-\mu _P} = {-\mu _Q}\) .
Given vector \( in \) d\(, the reflection of \)  under \((Q,P)\) ,
\( \text{reflect}(\mu _Q,\mu _P) = \mu _Q +(\mu _P -\) ,
satisfies:

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 [53]}, 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 [55]}, 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 [23]}, 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 [23]}, 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 [58]} 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 [59]}. 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 [3]} 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 [41]} 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 [3]} 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 [3]}.
To see why [3]} 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. [3]} 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 [3]} 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 [41]}), 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 [18]}. 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. [18]} 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 [18]}.
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 [53]}, 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 [72]}, 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 [73]}). 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 [74]}. 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 [74]} 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 [76]}, 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 \) 
[17]}
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 [18]} 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 [1]} - 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 [1]}, 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 [1]}.

Experiments
We run experiments on MNIST [82]}, 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>