# Machine Unlearning via Algorithmic Stability

Enayat Ullah Tung Mai Anup Rao Ryan Rossi Raman Arora Johns Hopkins
University. Email<EMAIL_ADDRESS>Research. Email<EMAIL_ADDRESS>Research. Email<EMAIL_ADDRESS>Research. Email<EMAIL_ADDRESS>Hopkins University. Email<EMAIL_ADDRESS>
###### Abstract

We study the problem of machine unlearning and identify a notion of
algorithmic stability, Total Variation (TV) stability, which we argue, is
suitable for the goal of exact unlearning. For convex risk minimization
problems, we design TV-stable algorithms based on noisy Stochastic Gradient
Descent (SGD). Our key contribution is the design of corresponding _efficient_
unlearning algorithms, which are based on constructing a (maximal) coupling of
Markov chains for the noisy SGD procedure. To understand the trade-offs
between accuracy and unlearning efficiency, we give upper and lower bounds on
excess empirical and populations risk of TV stable algorithms for convex risk
minimization. Our techniques generalize to arbitrary non-convex functions, and
our algorithms are differentially private as well.

###### Contents

  1. 1 Introduction
    1. 1.1 Our contributions
    2. 1.2 Related work
  2. 2 Problem setup and preliminaries
    1. 2.1 Streaming edit requests and exact unlearning
    2. 2.2 Convex risk minimization
    3. 2.3 Total variation stability and maximal coupling
  3. 3 Main results
  4. 4 Main ideas
    1. 4.1 Total variation stability from optimal transport
    2. 4.2 TV-stable learning algorithms and differential privacy
    3. 4.3 Unlearning via (un)couplings
  5. 5 Algorithms
    1. 5.1 TV-stable learning algorithms
    2. 5.2 Unlearning algorithms
  6. 6 Proofs of main results
    1. 6.1 Proof of Theorem 1
    2. 6.2 Proof of Theorem 2
    3. 6.3 Proof of Theorem 3
  7. 7 Discussion
  8. A Additional related work
  9. B Proofs for Section 5.1
  10. C Proofs for Section 5.2
    1. C.1 Unlearning for _sub-sample-GD_
    2. C.2 Unlearning for _noisy-m-A-SGD_
  11. D Runtime and space complexity
    1. D.1 Learning runtime
    2. D.2 Unlearning runtime
    3. D.3 Space complexity
  12. E Other algorithms and batch unlearning
    1. E.1 _noisy-m-SGD_
    2. E.2 _quantized-m-SGD_
  13. F Lower bounds on excess empirical risk
    1. F.1 Lower bound for mean computation
  14. G Excess population risk bounds
    1. G.1 Upper bounds
    2. G.2 Lower bounds
  15. H Algorithms for approximate unlearning
  16. I Experiments

##  1 Introduction

User data is employed in data analysis for various tasks such as drug
discovery, as well as in third-party services for tasks such as
recommendations. With such practices becoming ubiquitous, there is a growing
concern that sensitive personal data can get compromised. This has resulted in
a push for broader awareness of data privacy and ownership. These efforts have
led to several regulatory bodies enacting laws such as the European Union
General Data Protection Regulation (GDPR) and California Consumer Act, which
empowers the user with (among other things) the right to request to have their
personal data be _deleted_ (see Right to be forgotten, Wikipedia [31]).
However, currently it is unclear what it means to be _forgotten_ or have the
data _deleted_ in a rigorous sense. Nonetheless, there is a reasonable
expectation that merely deleting a user’s data from the database without
undoing the computations derived from the said data is insufficient. In the
settings where user data are directly utilized to build machine learning
models for, say prediction tasks, a reasonable criterion is that the system’s
state is adjusted to what it would have been if the user data were absent to
begin with – this is the criteria we adopt in our work. We refer to it as
_exact unlearning_ (see Definition 1 for a formal definition).

A straightforward way to comply with the requirement of exact unlearning is to
recompute (or retrain, in context of machine learning). This, however, is
often computationally expensive, and hence the goal is to design _efficient_
unlearning algorithms - herein and afterwards, we use efficient, in the
context of unlearning, to mean that the runtime is smaller than recompute
time. Most of the prior work focuses on specific structured problems (eg:
linear regression, clustering etc.) and the proposed unlearning algorithms
carefully leverage this structure for efficiency - see a discussion of useful
algorithmic principles like linearity, modularity etc. which enable efficient
unlearning, in [17]. In this work, we consider a larger class of problems:
smooth convex empirical risk minimization (ERM), which includes many machine
learning problems, eg: linear and logistic regression, as special cases (see
Section 2 for definitions). This class of smooth convex ERM is sufficiently
rich and arguably lacks structure useful in prior works. Gradient-based
optimization is a key algorithmic technique used for smooth convex ERM (and
most of machine learning and even beyond). Unlike prior works where the
(learning) algorithms for the specific problems (like linear regression) are
tailored enough to the problem to be amenable to efficient unlearning, an
optimization method hardly has any such useful structure. In particular, the
sequential nature of gradient descent makes it challenging to design non-
trivial unlearning algorithms, at least those which satisfy an _exact_
unlearning criterion. To elaborate, if the point to be deleted participates in
some iteration, then the subsequent steps are dependent on the to-be-deleted
point, and there is no known way but to redo the computations. It is natural
to then ask whether we can design unlearning algorithms with non-trivial
guarantees for this class of smooth convex ERM problems.

##### Problem statement (informal).

We consider a smooth convex ERM problem over a given initial dataset, in a
setting wherein we observe a stream of edits (insertion or deletion) to the
dataset. The goal is to design a learning algorithm that outputs an initial
model and a (corresponding) unlearning algorithm that updates the model after
an edit request. We require the following properties to hold: 1) exact
unlearning – at every time point in the stream, the output model is
indistinguishable from what we would have obtained if trained on the _updated_
dataset (i.e., without the deleted sample or with the inserted sample); 2) the
unlearning runtime should be _small_ ; 3) the output models should be
sufficiently accurate (measured in empirical risk).

###  1.1 Our contributions

##### Total variation stability.

We develop new algorithmic principles which enable exact unlearning in very
general settings. In particular, we propose a notion of algorithmic stability,
called total variation (TV) stability - an algorithmic property, which for any
problem, yields an _in-principle_ exact unlearning algorithm. Such an
algorithm might not be efficiently implementable computationally or due to the
data access restriction (sequential nature of edits). To demonstrate the
generality of our framework, we discuss, in Section E.2 how the previous work
of [17] for unlearning in $k$-means clustering using randomized quantization
can be interpreted as a special case of our framework - a TV stable method,
followed by efficient coupling based unlearning. We also note that the notion
of TV-stability has appeared before in [8], although in the seemingly
unrelated context of adaptive data analysis.

##### Convex risk minimization.

We make the above ideas of TV stability constructive in the special case of
smooth convex ERM problems. To elaborate, we give a TV stable learning
algorithm, and a corresponding efficient exact unlearning algorithm for smooth
convex ERM. Informally, for $n$ data points, and $d$ dimensional model and a
given $0<\rho\leq 1$, our method retrains only on $\rho$ fraction of edit
requests, while satisfying exact unlearning and maintaining that the accuracy
(excess empirical risk) is at most $\min\left\\{{\frac{1}{\sqrt{\rho
n}},\left({\frac{\sqrt{d}}{\rho n}}\right)^{4/5}}\right\\}$ (see Theorem 1 for
precise statement). This implies that for the (useful) regime of accuracy
greater than
$\min\left\\{{\frac{1}{\sqrt{n}},\left({\frac{\sqrt{d}}{n}}\right)^{4/5}}\right\\}$,
our algorithms provide a strict improvement over the only known method of re-
computation - see remarks after Theorem 1 for details. Furthermore, we also
give excess population risk bounds by leveraging known connections between
generalization and algorithmic stability (see Appendix G). Finally, we give
preliminary lower bounds on excess empirical and population risk for TV stable
algorithms for convex risk minimization.

##### Extensions.

Our results yield a number of interesting properties and extensions.

  * $\bullet$

Privacy: Even though privacy is not the goal of this work, some of our
$\rho$-TV stable algorithms, those based on noisy SGD like Algorithm 2, 5 are
$(\epsilon,\delta)$-differentially private (see Definition 3) with
$\epsilon=\rho\sqrt{\operatorname{log}\left(1/\delta\right)}$, for any
$\delta>0$. It is easy to see that these parameters can lie in the regime
reasonable for _good_ privacy properties i.e.
$\epsilon=O(1),\delta=\text{negl}(n)$. However, not all TV-stable algorithms,
for instance Algorithm 1, may have good privacy properties. Our work therefore
demonstrates interesting connections between techniques developed for
differential privacy and the problem of unlearning.

  * $\bullet$

Beyond Convexity: Interestingly, our unlearning techniques only require finite
sum structure in the optimization problem (for exact unlearning) and
Lipschitzness (for runtime bounds). Therefore our unlearning algorithms yield
provable unlearning for gradient-descent based methods for _any_ ERM or a
finite sum optimization problem. This means that we can apply the unlearning
algorithm even to non-convex problems, like training deep neural networks, and
everytime the unlearning algorithm does not recompute, it still guarantees
exact unlearning. As is typical, the accuracy in those cases is verified
empirically. Furthermore, Lipschitzness can be enforced by clipping gradients
- a popular heuristic in deep learning. In the worst-case (non-Lipschcitz)
scenario, it is easy to see we will need to recompute everytime, however this
is not the case in _typical_ situations - by which we mean when the deleted or
inserted data point is not an _outlier_ (measured in terms of its gradient
norm). This also means that our unlearning efficiency is based on, and can be
stated, in terms of, instance-dependent Lipschitz parameters, rather than the
worse case upper bound. In contrast, DP-training for non-convex models still
need Lipschitzness or clipping of gradients.

  * $\bullet$

Beyond gradient-based ERM: We also consider an approximate notion of
unlearning, based on differential privacy, which has appeared in the
literature [23, 16]. 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 Appendix H).

Our proofs are simple and conceptual and so we present the key ideas in the
main text and defer the proofs to the appendix.

###  1.2 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 [17] and [23], which
we discuss as follows. [17] 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 [17], 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 [17] as a special case (see
Section E.2 for details).

The work of [23] 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 [23] is approximate, based on
$(\epsilon,\delta)$-differential privacy, whereas we focus on exact
unlearning. Second, the unlearning runtime in [23] 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 [23] in Appendix A. We show therein that
with the same unlearning time, the accuracy guarantees of [23] are better than
us only in regimes where their approximate unlearning parameters and hence the
notion, is very weak.

##  2 Problem setup and preliminaries

###  2.1 Streaming edit requests and exact unlearning

We describe the setup very generally. Let ${\mathcal{Z}}$ be the data space,
$\Theta$ the output/parameter space, and ${\mathcal{M}}$ 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}:{\mathcal{Z}}^{*}\rightarrow\Theta\times{\mathcal{M}}$ is the
_batch_ learning algorithm, and
$\mathbf{U}:\Theta\times{\mathcal{M}}\times{\mathcal{Z}}\rightarrow\Theta\times{\mathcal{M}}$
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
${\mathcal{A}}(\cdot)$ denote the first output of $\mathbf{A}$ i.e.
${\mathcal{A}}(\cdot)=\mathbf{A}_{1}(\cdot)$. Similarly, let
${\mathcal{U}}(\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})=\left|{S\backslash
S^{\prime}}\right|+\left|{S^{\prime}\backslash S}\right|$. We now define exact
unlearning.

######  Definition 1 (Exact unlearning).

We say a procedure $(\mathbf{A},\mathbf{U})$ satisfies exact unlearning if for
any $S,S^{\prime}\subset{\mathcal{Z}}^{*}$ 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{\mathcal{M}}$, we have
$\mathbb{P}\left[{\mathbf{A}(S^{\prime})\in\mathcal{E}}\right]=\mathbb{P}\left[{\mathbf{U}(\mathbf{A}(S),S^{\prime}\backslash
S\cup S\backslash S^{\prime})\in\mathcal{E}}\right]$

######  Remark 1.

  1. 1.

A relaxation of the above definition is to maintain that only the output and
not the meta-data satisfy the above condition i.e.
${\mathcal{A}}(S^{\prime})={\mathcal{U}}(\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 [23].

  2. 2.

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}=\left\\{{\mathrm{z}_{1},\mathrm{z}_{2},\ldots,\mathrm{z}_{n}}\right\\},\mathrm{z}_{i}\in{\mathcal{Z}}$
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 [23], we
assume that at any point in the stream, the number of available data points is
at least $n/2$ and at most $2n$.

###  2.2 Convex risk minimization

We recall some basics from convex optimization. Let
${\mathcal{W}}\subset{\mathbb{R}}^{d}$ be a closed convex set such that
$\text{diameter}({\mathcal{W}})$ $\leq D$ where the diameter is measured in
Euclidean distance. Let ${\mathcal{Z}}$ be the instance space and let
$f:{\mathcal{W}}\times{\mathcal{Z}}\rightarrow{\mathbb{R}}$ be an
$L$-Lipschitz convex function in its first argument. For the constraint set
${\mathcal{W}}$, given a point $\mathrm{w}$, a projection function
${\mathcal{P}}:{\mathbb{R}}^{d}\rightarrow{\mathbb{R}}^{d}$ returns
${\mathcal{P}}(\mathrm{w})\in\arg\min_{\mathrm{v}\in{\mathcal{W}}}\left\|{\mathrm{w}-\mathrm{v}}\right\|$.
The function $f$ is $L$-smooth in its first argument if
$\left\|{\nabla_{\mathrm{w}}f(\mathrm{w}_{1},\mathrm{z})-\nabla_{\mathrm{w}}f(\mathrm{w}_{2},\mathrm{z})}\right\|\leq
L\left\|{\mathrm{w}_{1}-\mathrm{w}_{2}}\right\|\
\forall\mathrm{w}_{1},\mathrm{w}_{2}\in{\mathcal{W}},\mathrm{z}\in{\mathcal{Z}}$.
We will drop the subscript $\mathrm{w}$ 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=\left\\{{\mathrm{z}_{1},\mathrm{z}_{2},\ldots,\mathrm{z}_{n}}\right\\}$, we
look at the following class of problems, known as empirical risk minimization
(ERM).

$\displaystyle\min_{\mathrm{w}\in{\mathcal{W}}}\left\\{{\widehat{F}_{S}(\mathrm{w}):=\frac{1}{n}\sum_{j=1}^{n}f(\mathrm{w},\mathrm{z}_{i})}\right\\}$
(1)

Let ${\mathcal{A}}(S)$ be the output of algorithm ${\mathcal{A}}$ on dataset
$S$. We will give guarantees on _expected excess empirical risk_ , which is
$\mathbb{E}{\widehat{F}({\mathcal{A}}(S))}-F(\mathrm{w}_{S}^{*})$, where
$\mathrm{w}^{*}_{S}$ is the minimizer:
$\mathrm{w}^{*}_{S}\in\arg\min_{\mathrm{w}\in{\mathcal{W}}}\widehat{F}_{S}(\mathrm{w})$
and the expectation is taken with respect to the randomness in algorithm
${\mathcal{A}}$. 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 Appendix G.

###  2.3 Total variation stability and maximal coupling

We first state the definition of Total Variation (TV) distance between two
distributions $P$ and $Q$.

$\displaystyle\text{TV}(P,Q)=\sup_{\text{measurable sets
}R}\left|{P(R)-Q(R)}\right|=\frac{1}{2}\left\|{\phi_{P}-\phi_{Q}}\right\|_{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.

######  Definition 2 ($\rho$-TV-stability).

An algorithm ${\mathcal{A}}$ is said to be $\rho$-TV-stable if

$\displaystyle\sup_{S,S^{\prime}:\Delta(S,S^{\prime})=1}\text{TV}({\mathcal{A}}(S),{\mathcal{A}}(S^{\prime}))\leq\rho$

######  Remark 2.

  1. 1.

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

  2. 2.

Suppose $S$ is a dataset of $n$ points, and $S^{\prime}$ is a dataset of
$n+k_{2}$ points such that $\left|{S\backslash S^{\prime}}\right|=k_{1}$.
Then, if algorithm ${\mathcal{A}}$ is $\rho$-TV stable, then by triangle
inequality of TV and repeated applications of the above definition, we have
that
$\text{TV}({\mathcal{A}}(S),{\mathcal{A}}(S^{\prime}))\leq(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 $({\mathcal{X}},{\mathcal{B}})$, where ${\mathcal{B}}$
denotes the sigma-algebra on ${\mathcal{X}}$, is a distribution
$\pi\in{\mathbb{P}}({\mathcal{X}}\times{\mathcal{X}},{\mathcal{B}}\otimes{\mathcal{B}})$
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.

  1. 1.

For any coupling $\pi\in\Pi(P,Q)$, if the random variable $(p,q)\sim\pi$, then
$\text{TV}(P,Q)\leq\mathbb{P}\left[{p\neq q}\right]$.

  2. 2.

There exists a “maximal" coupling $\pi^{*}$ such that if $(p,q)\sim\pi^{*}$,
then $\text{TV}(P,Q)=\mathbb{P}\left[{p\neq q}\right]$

The above establishes that
$\text{TV}(P,Q)=\inf_{\pi\in\Pi(P,Q)}\mathbb{P}_{(p,q)\sim\pi}\left[p\neq
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.

##  3 Main results

We state our main result on designing learning and unlearning algorithms in a
stream of edit requests.

######  Theorem 1 (Main Theorem).

For any $\frac{1}{n}\leq\rho<\infty$, there exist a learning and a
corresponding unlearning algorithm such that for any $f(\cdot,\mathrm{z})$,
which is $L$-smooth and $G$-Lipschitz convex function $\forall\ \mathrm{z}$,
and a stream of edit requests,

  1. 1.

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

  2. 2.

At time $i$ in the stream, outputs $\widehat{\mathrm{w}}_{S^{i}}$ with excess
empirical risk bounded as,

$\displaystyle\mathbb{E}{\widehat{F}_{S^{i}}(\widehat{\mathrm{w}}_{S^{i}})-\widehat{F}_{S^{i}}(\mathrm{w}_{S^{i}}^{*})}\lesssim\min\left\\{{\frac{GD}{\sqrt{\rho
n}},\left({\frac{L^{1/4}GD^{3/2}\sqrt{d}}{(\rho n)}}\right)^{4/5}}\right\\}$

  3. 3.

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

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 Section D.1).

##### 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\left\\{{\frac{GD}{\sqrt{n}},\left({\frac{L^{1/4}GD^{3/2}\sqrt{d}}{n}}\right)^{4/5}}\right\\}$.
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 Appendix G
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 Theorem 1 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 Section 5.

##### 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 Section 4.1 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
$\left({\frac{\sqrt{d}}{\rho n}}\right)^{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 $\left({\frac{\sqrt{d}}{\rho n}}\right)^{2/3}$ (see
Section E.1 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 Theorem 1 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.

######  Theorem 2 (Upper bound).

For any $0<\rho<\infty$, there exists an algorithm which is
$\min\left\\{{\rho,1}\right\\}$-TV stable, such that for any
$f(\cdot,\mathrm{z})$ which is $L$-smooth and $G$-Lipschitz convex function
$\forall\ \mathrm{z}$, and any dataset $S$ of $n$ points, outputs
$\widehat{\mathrm{w}}_{S}$ which satisfies the following.

$\displaystyle\mathbb{E}{\widehat{F}_{S}(\widehat{\mathrm{w}}_{S})-\widehat{F}_{S}(\mathrm{w}^{*}_{S})}\lesssim
GD\min\left\\{{\frac{1}{\sqrt{\rho n}},\frac{\sqrt{d}}{\rho n}}\right\\}$

We show that the condition $\rho\geq\frac{1}{n}$ in Theorem 2 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\geq 1$ in
our lower bound since it puts no constraint on the algorithm.

######  Theorem 3 (Lower bound).

For any $\rho$-TV-stable algorithm ${\mathcal{A}}$, 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:

  1. 1.

For any $0<\rho<1$, and any dimension $d$,
$\mathbb{E}{\widehat{F}_{S}({\mathcal{A}}(S))-\widehat{F}_{S}(\mathrm{w}_{S}^{*})}\gtrsim
GD\min\left\\{{1,\frac{1}{\rho n}}\right\\}$

  2. 2.

Assuming that ${\mathcal{A}}(S)$ has a probability density function upper
bounded by $K\leq O(2^{d})$, then for
$n>72,\frac{1}{n}\leq\rho\leq\frac{1}{4}$ and large enough $d$,
$\mathbb{E}{\widehat{F}_{S}({\mathcal{A}}(S))-\widehat{F}_{S}(\mathrm{w}_{S}^{*})}\gtrsim
GD\min{\left\\{{1,\frac{1}{\sqrt{\rho n}}}\right\\}}$

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 Theorem 2 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}\leq\frac{1}{\sqrt{\rho n}}\iff d\leq\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 ${\mathcal{A}}(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 Appendix F for derivations of the lower bounds.

##  4 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.

###  4.1 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={\mathcal{A}}(S)$ and $Q={\mathcal{A}}(S^{\prime})$
for some randomized algorithm ${\mathcal{A}}$. 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
[28], which briefly explain below. Given probability distributions $P$ and $Q$
over measurable space ${\mathcal{X}}$, and a cost function
$c:{\mathcal{X}}\times{\mathcal{X}}\rightarrow{\mathbb{R}}$, 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)=\begin{cases}1&\text{ if }x\neq y\\\ 0&\text{
otherwise}\end{cases}$. 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}\left\\{{x\neq
y}\right\\}$ \- 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.

###  4.2 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 [12], defined as follows.

######  Definition 3 (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
${\mathcal{E}}\in\text{Range}({\mathcal{A}})$,

$\displaystyle\mathbb{P}(\mathcal{A}(S)\in{\mathcal{E}})\leq
e^{\epsilon}\mathbb{P}(\mathcal{A}(S^{\prime})\in{\mathcal{E}})+\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{\operatorname{log}\left(1/\delta\right)}}{\epsilon}$ to
a $1$-sensitive function, provides $(\epsilon,\delta)$-DP [14]. It can be
shown that the same method also provides $\rho$-TV stability, with
$\rho=\frac{\epsilon}{\sqrt{\operatorname{log}\left(1/\delta\right)}}$.

##### 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 [9] 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 $\mathrm{w}_{j}$ denoted as
$\nabla\widehat{F}_{S}(\mathrm{w}_{j},\mathrm{z}_{b_{j}})$ and update as
follows:

$\mathrm{w}_{j+1}=\mathrm{w}_{j}-\eta\left({\nabla\widehat{F}_{S}(\mathrm{w}_{j},\mathrm{z}_{b_{j}})+\theta_{t}}\right)$

where $\theta_{j}\sim{\mathcal{N}}(0,\sigma^{2}{\mathbb{I}}_{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 1.

Mini-
batches$b_{1}$$b_{2}$$b_{3}$$b_{T}$$\mathrm{w}_{2}$$\mathrm{w}_{1}$$\mathrm{w}_{3}$$\mathrm{w}_{4}$$\mathrm{w}_{T+1}$Iterates……
Figure 1: Markov chain for the _noisy-m-SGD_ algorithm

##### 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 [23], 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 [9], 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 Section 4.3.

###  4.3 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={\mathcal{A}}(S)$ to $Q={\mathcal{A}}(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 $\left\\{{(p,q):p=q}\right\\}$ where
$(p,q)\sim\pi$, and similarly, the non-diagonal is the set
$\left\\{{(p,q):p\neq q}\right\\}$, $(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}\left\\{{(p,q):p\neq
q}\right\\}=\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 1. 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}\left({\left\\{{\mathrm{w}_{j}^{P}}\right\\}_{j=1}^{T+1},\left\\{{\mathrm{w}_{j}^{Q}}\right\\}_{j=1}^{T+1}}\right)\leq\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 [29] and [15], 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 1. 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 $\left\\{{b^{P}_{j}}\right\\}_{j=1}^{T}$ and
$\left\\{{b^{Q}_{j}}\right\\}_{j=1}^{T}$. We then compute an estimate of
marginal densities of $\mathrm{w}_{j}$ under $P$ and $Q$, via conditional
densities under coupled mini-batches i.e. we compute
$\phi_{P}(\mathrm{w}_{j}|b_{j}^{Q})$ and $\phi_{Q}(\mathrm{w}_{j}|b_{j}^{P})$
\- note that these are just Gaussian densities evaluated at
$\mathrm{w}_{j},b_{j}^{P}$ and $\mathrm{w}_{j},b_{j}^{Q}$ respectively, where
$\mathrm{w}_{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\leq\frac{\phi_{Q}(\mathrm{w}_{j}|b_{j}^{Q})}{\phi_{P}(\mathrm{w}_{j}|b_{j}^{P})}$.
If the step succeeds, we accept $\mathrm{w}_{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 $\mathrm{w}_{t+1}^{Q}$ by _reflecting_
$\mathrm{w}_{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 Algorithm 4
(please see Section 5.2.2 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}(\mathrm{w}_{j}|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 Theorem 3). Please see Section D.2 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 Section D.1, we also remark that
using even faster algorithms like Katyusha [2] 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 $\mathrm{w}_{j}$ and $\mathrm{w}_{j}^{\prime}$. The distribution
of $\mathrm{w}_{j}$ is such that, for any event $E$ in range of
$\mathrm{w}_{j}$, $P(\mathrm{w}_{j}\in E)=\mathbb{E}_{b}P(\mathrm{w}_{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}(\mathrm{w}_{j})}{\phi_{P}(\mathrm{w}_{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}(\mathrm{w}_{j}|b_{j}^{Q})}{\phi_{P}(\mathrm{w}_{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}((\mathrm{w}^{P}_{j},b^{P}_{j}),(\mathrm{w}^{Q}_{j},b^{Q}_{j}))\geq\mathbb{E}_{(b_{j}^{P},b_{j}^{Q})}\text{TV}(\mathrm{w}_{j}^{P}|b_{j}^{P},\mathrm{w}_{j}^{Q}|b_{j}^{Q})\geq\text{TV}(\mathrm{w}_{j}^{P},\mathrm{w}_{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.

##  5 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 Appendix D, we
explain what data structures to use for computational efficiency. The proofs
of results in this section are deferred to Appendix B.

###  5.1 TV-stable learning algorithms

####  5.1.1 _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 (Algorithm 3) uses all the saved iterates.
However this is done only for ease of presentation - in Section D.3, we
discuss a simple efficient implementation (of the unlearning algorithm), which
doesn’t need _any_ iterate, yet has the same unlearning time complexity.

Algorithm 1 _sub-sample-GD_ $(\mathrm{w}_{t_{0}},t_{0})$

0: Initial model $\mathrm{w}_{t_{0}}$, data points
$\left\\{{\mathrm{z}_{1},\ldots,\mathrm{z}_{n}}\right\\},T,m,\eta$

1: for $t=t_{0},t_{0}+1\ldots,T$ do

2: Sample mini-batch $b_{t}$ of size $m$ uniformly randomly

3: $\mathrm{g}_{t}=\frac{1}{m}\sum_{j\in b_{t}}\nabla
f(\mathrm{w}_{t},\mathrm{z}_{j})$

4:
$\mathrm{w}_{t+1}={\mathcal{P}}\left({\mathrm{w}_{t}-\eta\mathrm{g}_{t}}\right)$

5: Save($b_{t},\mathrm{w}_{t},\mathrm{g}_{t}$)

6: end for

6: $\widehat{\mathrm{w}}_{S}=\frac{1}{T}\sum_{t=1}^{T+1}\mathrm{w}_{t}$

We now give guarantees on excess empirical risk for _sub-sample-GD_.

######  Proposition 1.

Let $f(.,\mathrm{z})$ be an $L$-smooth $G$-Lipschitz convex function $\forall\
\mathrm{z}$. Algorithm 1, run with
$t_{0}=1,\eta=\min\left\\{{\frac{1}{2L},\frac{D\sqrt{\rho
n}}{GT}}\right\\},T=\frac{DL\sqrt{\rho n}}{G}$, and
$m=\max{\left\\{{\frac{G\sqrt{\rho n}}{DL},1}\right\\}}$, outputs
$\widehat{\mathrm{w}}_{S}$ which is $\min\left\\{{\rho,1}\right\\}$-TV-stable
and satisfies
$\mathbb{E}{\widehat{F}_{S}(\widehat{\mathrm{w}}_{S})-\widehat{F}_{S}(\mathrm{w}^{*}_{S})}\lesssim\frac{GD}{\sqrt{\rho
n}}.$

####  5.1.2 _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.

Algorithm 2 _noisy-m-A-SGD_($\mathrm{w}_{t_{0}},t_{0})$

0: Initial model $\mathrm{w}_{t_{0}}$, data points
$\left\\{{\mathrm{z}_{1},\ldots,\mathrm{z}_{n}}\right\\}$, $T$, $\eta$, $m$

1: $\mathrm{w}_{0}=0$

2: for $t=t_{0},t_{0}+1\ldots,T$ do

3: Sample mini-batch $b_{t}$ of size $m$ uniformly randomly

4: Sample $\theta_{t}\sim{\mathcal{N}}(0,\sigma^{2}{\mathbb{I}}_{d})$

5:
$\mathring{\mathrm{w}}_{t}=(1-\alpha_{t})\mathrm{w}_{t}+\alpha_{t}\mathrm{w}_{t-1}$

6: $\mathrm{g}_{t}=\frac{1}{m}\sum_{j\in b_{t}}\nabla
f(\mathring{\mathrm{w}}_{t},\mathrm{z}_{j})$

7:
$\mathrm{w}_{t+1}={\mathcal{P}}\left({\mathring{\mathrm{w}}_{t}-\eta\left({\mathrm{g}_{t}+\theta_{t}}\right)}\right)$

8:
Save($b_{t},\theta_{t},\mathrm{w}_{t},\mathring{\mathrm{w}}_{t},\mathrm{g}_{t}$)

9: end for

9: $\widehat{\mathrm{w}}_{S}=\mathrm{w}_{T+1}$

We now state our results for _noisy-m-A-SGD_.

######  Proposition 2.

Let $f(.,\mathrm{z})$ be an $L$-smooth $G$-Lipschitz convex function $\forall\
\mathrm{z}$. For any $0<\rho<\infty$, Algorithm 2, run with
$t_{0}=1,\eta=\min\left\\{{\frac{1}{2L},\frac{D}{\left({\frac{G}{\sqrt{m}}+\sigma}\right)T^{3/2}}}\right\\}$,
$\alpha_{0}=0,\alpha_{t}=\frac{1-t}{t+2}$, $\sigma=\frac{8\sqrt{T}G}{n\rho}$,
and $T\geq\frac{(n\rho)^{2}}{16m^{2}}$ outputs $\widehat{\mathrm{w}}_{S}$
which is $\min\left\\{{\rho,1}\right\\}$-TV stable and satisfies

$\mathbb{E}{\widehat{F}_{S}(\widehat{\mathrm{w}}_{S})-\widehat{F}_{S}(\mathrm{w}^{*}_{S})}\lesssim\frac{LD^{2}}{T^{2}}+\frac{GD}{\sqrt{Tm}}+\frac{GD\sqrt{d}}{n\rho}.$

Choosing $T$ and $m$ appropriately gives us the following corollary.

######  Corollary 1.

Let $f(.,\mathrm{z})$ be an $L$-smooth $G$-Lipschitz convex function $\forall\
\mathrm{z}$. For any $0<\rho<\infty$, Algorithm 2, run with
$t_{0}=1,m\geq\min\left\\{{\frac{d}{16},\frac{1}{4}\left({\frac{(\rho
n)^{3}G\sqrt{d}}{LD}}\right)^{1/4}}\right\\}$,
$\eta=\min\left\\{{\frac{1}{2L},\frac{D}{\left({\frac{G}{\sqrt{m}}+\sigma}\right)T^{3/2}}}\right\\}$,
$\alpha_{0}=0,\alpha_{t}=\frac{1-t}{t+2}$, $\sigma=\frac{8\sqrt{T}G}{n\rho}$,
and $T=\max\left\\{{\frac{(\rho n)^{2}}{md},\sqrt{\frac{LD\rho
n}{G\sqrt{d}}}}\right\\}$ outputs $\widehat{\mathrm{w}}_{S}$ which is
$\min\left\\{{\rho,1}\right\\}$-TV stable and satisfies
$\mathbb{E}{\widehat{F}_{S}(\widehat{\mathrm{w}})-\widehat{F}_{S}(\mathrm{w}^{*}_{S})}\lesssim\frac{GD\sqrt{d}}{\rho
n}$.

######  Remark 3.

The choice of $T$ in Corollary 1 yields that the largest mini-batch size that
can be set, without hurting runtime, is $m=\left({\frac{(\rho
n)^{3}G}{\sqrt{d}^{3}LD}}\right)^{1/2}=\left({\frac{G}{LD}}\right)^{2}T^{3}$.
Furthermore, the condition
$m\geq\min\left\\{{\frac{d}{16},\frac{1}{4}\left({\frac{(\rho
n)^{3}G\sqrt{d}}{LD}}\right)^{1/4}}\right\\}$ yields $(\rho
n)\geq\left({\frac{LD(\sqrt{d})^{7}}{256G}}\right)^{1/3}$.

In Proposition 5 (in Appendix , Section 5.1), we show that the upper bound on
total variation stability parameter of Algorithm 2 derived in Proposition 2 is
tight in all problem parameters, upto constants.

###  5.2 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 1
and the other for Algorithm 2, and show that the probability with which a
recompute is triggered is small - please see Appendix D, for a finer analysis
of runtime. The proofs of results in this section are deferred to Appendix C.

####  5.2.1 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 Algorithm 3 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 3 implements the above
procedure.

Algorithm 3 Unlearning for _sub-sample-GD_

0: Data point index $j$ to delete or data point $\mathrm{z}$ to insert (index
$n+1$)

1: for $t=1,2\ldots,T$ do

2: Load$\left({b_{t},\mathrm{g}_{t},\mathrm{w}_{t}}\right)$

3: if deletion and $j\in b_{t}$ then

4:  _sub-sample-GD_ $(\mathrm{w}_{t},t)$ // Continue training on current
dataset

5: break

6: else if insertion and Bernoulli$\left({\frac{m}{n+1}}\right)$ then

7: Sample $i\sim\text{Uniform}(b_{t})$

8: $\mathrm{g}_{t}^{\prime}=\mathrm{g}_{t}-\frac{1}{m}\left({\nabla
f(\mathrm{w}_{t},\mathrm{z}_{i})-\nabla f(\mathrm{w}_{t},\mathrm{z})}\right)$

9:
$\mathrm{w}_{t+1}={\mathcal{P}}\left({\mathrm{w}_{t}-\eta\left({\mathrm{g}_{t}^{\prime}+\theta_{t}}\right)}\right)$

10:
Save($\mathrm{w}_{t+1},\mathrm{g}_{t}^{\prime},b_{t}\backslash\left\\{{i}\right\\}\cup\left\\{{n+1}\right\\}$)

11:  _sub-sample-GD_ $(\mathrm{w}_{t+1},t+1)$ // Continue training on current
dataset

12: break

13: end if

14: end for

We state our main result for unlearning with Algorithm 3 below.

######  Proposition 3.

(Algorithm 1, Algorithm 3) satisfies exact unlearning. Moreover, for $k$
edits, Algorithm 3 recomputes with probability at most $2k\rho$.

####  5.2.2 Unlearning for _noisy-m-A-SGD_

Algorithm 4 Unlearning for _noisy-m-A-SGD_

0: Data point index $j$ to delete or data point $\mathrm{z}$ to insert (index
$n+1$)

1: for $t=1,2\ldots,T$ do

2:
Load$\left({\theta_{t},\mathrm{w}_{t},\mathring{\mathrm{w}}_{t},b_{t},\mathrm{g}_{t}}\right)$

3: if deletion and $j\in b_{t}$ then

4: Sample $i\sim\text{Uniform}([n]\backslash b_{t})$

5: $\mathrm{g}_{t}^{\prime}=\mathrm{g}_{t}-\frac{1}{m}\left({\nabla
f(\mathring{\mathrm{w}}_{t},\mathrm{z}_{j})-\nabla
f(\mathring{\mathrm{w}}_{t},\mathrm{z}_{i})}\right)$

6:
Save($\mathrm{g}_{t}^{\prime},b_{t}\backslash\left\\{{j}\right\\}\cup\left\\{{i}\right\\}$)

7: else if insertion and Bernoulli$\left({\frac{m}{n+1}}\right)$ then

8: Sample $i\sim\text{Uniform}(b_{t})$

9: $\mathrm{g}_{t}^{\prime}=\mathrm{g}_{t}-\frac{1}{m}\left({\nabla
f(\mathring{\mathrm{w}}_{t},\mathrm{z}_{i})-\nabla
f(\mathring{\mathrm{w}}_{t},\mathrm{z})}\right)$

10:
Save($\mathrm{g}_{t}^{\prime},b_{t}\backslash\left\\{{i}\right\\}\cup\left\\{{n+1}\right\\}$)

11: else

12: continue

13: end if

14: $\xi_{t}=\mathrm{g}_{t}+\theta_{t}$

15: if
Uniform$\left({0,1}\right)\geq\frac{\phi_{{\mathcal{N}}(\mathrm{g}_{t}^{\prime},\sigma^{2}{\mathbb{I}})}(\xi_{t})}{\phi_{{\mathcal{N}}(\mathrm{g}_{t},\sigma^{2}{\mathbb{I}})}(\xi_{t})}$
then

16:
$\xi^{\prime}_{t}=\text{reflect}(\xi_{t},\mathrm{g}_{t}^{\prime},\mathrm{g}_{t})$

17: $\mathrm{w}_{t+1}=\mathrm{w}_{t}-\eta\xi_{t}^{\prime}$

18: $\text{Save}(\xi^{\prime}_{t})$

19:  _noisy-m-A-SGD( $\mathrm{w}_{t+1},t+1$)_ // Continue retraining on
current dataset

20: break

21: end if

22: end for

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
Section 4.3, and we just fill in some details here. We first describe how
Algorithm 4 couples mini-batch indices while handling edit request.

##### Coupling mini-batch indices:

After observing a deletion request, in Algorithm 4, 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.

######  Definition 4 (Reflection map).

Given a vector $\mathrm{u}$ and two vectors $\mathrm{x}$ and $\mathrm{y}$, the
reflection of $\mathrm{u}$ under $(\mathrm{x},\mathrm{y})$, denoted as
$\text{reflect}(\mathrm{u},\mathrm{x},\mathrm{y})$, is defined as

$\displaystyle\text{reflect}(\mathrm{u},\mathrm{x},\mathrm{y})=\mathrm{x}+(\mathrm{y}-\mathrm{u})$

Reflection coupling is a classical idea in probability, used to construct
couplings between symmetric probability distributions [21]. The reflection
map, given $\mathrm{u},\mathrm{x},\mathrm{y}$, reflects $\mathrm{u}$ about the
mid-point of $\mathrm{x}$ and $\mathrm{y}$. The context in which we will use
it is $\mathrm{u}$ would be a sampled point from a Gaussian under old dataset
$S$ (on which the model was trained on), and $\mathrm{x}$ and $\mathrm{y}$
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 Section C.2.2 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
$\mathrm{w}_{t+1}$ using rejection sampling. To elaborate, at each iteration,
we check if the noisy iterate, defined as
$\bar{\mathrm{w}}_{t+1}=\mathring{\mathrm{w}}_{t}-\eta(\mathrm{g}_{t}+\theta_{t})$
is a _good_ sample for the dataset $S^{\prime}$, where $\mathrm{g}_{t}$ is the
gradient computed on $\mathring{\mathrm{w}}_{t}$ using a uniform sub-sample
from $S$. To do this, we need to compute a ratio of _estimated_ marginal
densities of $\mathrm{w}_{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{\mathrm{w}}_{t+1}=\mathring{\mathrm{w}}_{t}-\eta(\mathrm{g}_{t}+\theta_{t})$.
If we also fix the sampled mini-batch $b_{t}$, then $\bar{\mathrm{w}}_{t+1}$
is distributed as
${\mathcal{N}}(\mathring{\mathrm{w}}_{t}-\eta\mathrm{g}_{t},\eta^{2}\sigma^{2}{\mathbb{I}})$.
However, once we unfix $b_{t}$, then $\mathrm{w}_{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 $\mathrm{w}_{t+1}$ (and $\mathrm{w}_{t+1}^{\prime}$ \- the
iterate for dataset $S^{\prime}$) evaluated at
$\mathring{\mathrm{w}}_{t}-\eta(\mathrm{g}_{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 Algorithm 2, 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 Section C.2.3 for a more formal treatment of the coupling
procedure. We now state the main result for this section.

######  Proposition 4.

(Algorithm 2, Algorithm 4) satisfies exact unlearning. Moreover, for $k$
edits, Algorithm 4 recomputes with probability at most
$\frac{k\rho\sqrt{T}}{4}$

##  6 Proofs of main results

In this section, we give the proofs of main results, stated in Section 3,
using the results in the preceding sections.

###  6.1 Proof of Theorem 1

The proof follows by combining the guarantees for the two algorithms we
present: _sub-sample-GD_ (Algorithm 1) and _noisy-m-A-SGD_ (Algorithm 2), and
their corresponding unlearning algorithms: Algorithm 3 and Algorithm 4. We
discuss these one by one. From Proposition 1, we have that, given $0<\rho\leq
1$, _sub-sample-GD_ is $\rho$-TV stable and has excess empirical risk bounded
by $O\left({\frac{GD}{\sqrt{\rho n}}}\right)$. This holds at every point in
the stream by assumption that the number of samples are between $\frac{n}{2}$
and $2n$. Furthermore, from Proposition 3, 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 7,
we have that there exist efficient implementations, such that the runtime of
unlearning for _sub-sample-GD_ is
$O(\max\left\\{{k,\min\left\\{{\rho,1}\right\\}k\cdot\text{Training
time}}\right\\}$, 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
Proposition 2, for dataset $S$ of $n$ points, we have that, given
$0<\tilde{\rho}\leq 1$, _noisy-m-A-SGD_ is $\tilde{\rho}$-TV stable and its
excess empirical risk is bounded as follows:

$\mathbb{E}{\widehat{F}_{S}(\widehat{\mathrm{w}}_{S})-\widehat{F}_{S}(\mathrm{w}^{*}_{S})}\lesssim\frac{LD^{2}}{T^{2}}+\frac{GD}{\sqrt{Tm}}+\frac{GD\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 4, 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
8, we have that there exist efficient implementations, such that the runtime
of unlearning for _noisy-m-A-SGD_ is
$O(\max\\{k,k\min\left\\{{\tilde{\rho}\sqrt{T},1}\right\\}\cdot\text{Training
time}\\})$. In the statement of Theorem 1, 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:

$\mathbb{E}{\widehat{F}_{S}(\widehat{\mathrm{w}}_{S})-\widehat{F}_{S}(\mathrm{w}^{*}_{S})}\lesssim\frac{LD^{2}}{T^{2}}+\frac{GD}{\sqrt{Tm}}+\frac{GD\sqrt{d}\sqrt{T}}{n\rho}.$

We use the largest mini-batch size, which does not hurt runtime, which is
$m=\left({\frac{G}{LD}}\right)^{2}T^{3}$. This simplifies the upper bound to
$\frac{LD^{2}}{T^{2}}+\frac{GD\sqrt{d}\sqrt{T}}{n\rho}$. Optimizing the trade-
off, we have $\frac{LD^{2}}{T^{2}}=\frac{GD\sqrt{d}\sqrt{T}}{n\rho}\iff
T=\left({\frac{LD(n\rho)}{G\sqrt{d}}}\right)^{2/5}$, and the excess empirical
risk becomes
$\mathbb{E}{\widehat{F}_{S}(\widehat{\mathrm{w}}_{S})-\widehat{F}_{S}(\mathrm{w}^{*}_{S})}\lesssim\frac{LD^{2}}{T^{2}}=\left({\frac{L^{1/4}GD^{3/2}\sqrt{d}}{(\rho
n)}}\right)^{4/5}$ – this recovers the other term in the upper bound in
Theorem 1. However, note that Proposition 2 has an additional condition that
$T\geq\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=\left({\frac{G}{LD}}\right)^{2}T^{3}$ and
$\tilde{\rho}=\rho/\sqrt{T}$, the condition simplifies as
$T\geq\frac{(n\rho)^{2}}{16T\cdot T^{6}}\left({\frac{LD}{G}}\right)^{2}\iff
T^{8}\geq\frac{(n\rho)^{2}}{16}\left({\frac{LD}{G}}\right)^{4}\iff\left({\frac{LD(n\rho)}{G\sqrt{d}}}\right)^{16/5}\geq\frac{(n\rho)^{2}}{16}\left({\frac{LD}{G}}\right)^{4}\iff\left({\frac{\sqrt{d}}{(n\rho)}}\right)^{4/5}\leq\frac{2}{\sqrt{n\rho}}\left({\frac{LD}{G}}\right)^{1/5}\iff\left({\frac{L^{1/4}GD^{3/2}\sqrt{d}}{(\rho
n)}}\right)^{4/5}\leq\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 Theorem 1. ∎

###  6.2 Proof of Theorem 2

We give two algorithms, _sub-sample-GD_ (Algorithm 1) and _noisy-m-A-SGD_
(Algorithm 2), one for each of the upper bounds. From Proposition 1 and
Corollary 1, we have that, given $0<\rho<\infty$, these are
$\min\left\\{{\rho,1}\right\\}$-TV stable and their excess empirical risk is
bounded is $O\left({\frac{GD}{\sqrt{\rho n}}}\right)$ and
$O\left({\frac{GD\sqrt{d}}{\tilde{\rho}n}}\right)$ respectively. Hence
combining the above by taking a minimum, establishes the claimed result. ∎

###  6.3 Proof of Theorem 3

In all the lower bounds, we have a $GD$ term - this is a trivial lower bound,
since if an algorithm is defined as ${\mathcal{A}}(S)=0$ (or any constant),
then this is perfectly $TV$ stable ($\rho=0$), and the expected excess
empirical risk is upper bounded as
$\widehat{F}_{S}({\mathcal{A}}(S))-\widehat{F}_{S}(\mathrm{w}_{S}^{*})\leq
G\left\|{{\mathcal{A}}(S)-\mathrm{w}_{S}^{*}}\right\|\leq GD$, where the first
inequality uses $G$-Lipschitzness of $\widehat{F}_{S}$ and the second the fact
the both ${\mathcal{A}}(S)$ and $\mathrm{w}_{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 [9], we consider $G=D=1$, since a simple reduction
gives a factor of $GD$ for general $G$ and $D$. Furthermore, similar to [9],
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 Proposition 10. 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}\leq\left\|{\mu(S)}\right\|\leq 2M$. The accuracy, denoted by
$\alpha$, is defined as
$\alpha^{2}=\mathbb{E}\left\|{{\mathcal{A}}(S)-\mu(S)}\right\|^{2}$,
${\mathcal{A}}$ is a $\rho$-TV stable algorithm, and the expectation is taken
over the algorithm’s randomness. The first part of Theorem 3 follows Theorem 6
which is based on a simple reduction argument. This gives us that
$\alpha\geq\frac{1}{\rho n}$ with $M=\frac{1}{\rho n}$. Plugging it in
Proposition 10, this gives us that excess empirical risk is lower bounded by
$\Omega\left({\frac{1}{\rho n}}\right)$. Similarly, the second part follows
from Theorem 7 which gives us $\alpha\geq\frac{1}{\sqrt{\rho n}}$ with
$M=\frac{1}{\sqrt{\rho n}}$ \- the condition $\alpha\leq\frac{1}{4}$ in the
statement of Theorem 7 can be absorbed in the trivial lower bound $GD$. ∎

##  7 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.

  1. 1.

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.

  2. 2.

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.

  3. 3.

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.

## References

  * ACG+ [16] Martin Abadi, Andy Chu, Ian Goodfellow, H Brendan McMahan, Ilya Mironov, Kunal Talwar, and Li Zhang.  Deep learning with differential privacy.  In Proceedings of the 2016 ACM SIGSAC Conference on Computer and Communications Security, pages 308–318, 2016.
  * AZ [17] Zeyuan Allen-Zhu.  Katyusha: The first direct acceleration of stochastic gradient methods.  The Journal of Machine Learning Research, 18(1):8194–8244, 2017\.
  * AZ [18] Zeyuan Allen-Zhu.  How to make the gradients small stochastically: Even faster convex and nonconvex sgd.  In Advances in Neural Information Processing Systems, pages 1157–1167, 2018.
  * BBG [18] Borja Balle, Gilles Barthe, and Marco Gaboardi.  Privacy amplification by subsampling: Tight analyses via couplings and divergences.  In Advances in Neural Information Processing Systems, pages 6277–6287, 2018.
  * BCCC+ [19] Lucas Bourtoule, Varun Chandrasekaran, Christopher Choquette-Choo, Hengrui Jia, Adelin Travers, Baiwu Zhang, David Lie, and Nicolas Papernot.  Machine unlearning.  arXiv preprint arXiv:1912.03817, 2019.
  * BE [02] Olivier Bousquet and André Elisseeff.  Stability and generalization.  Journal of machine learning research, 2(Mar):499–526, 2002.
  * BFTT [19] Raef Bassily, Vitaly Feldman, Kunal Talwar, and Abhradeep Guha Thakurta.  Private stochastic convex optimization with optimal rates.  In Advances in Neural Information Processing Systems, pages 11282–11291, 2019.
  * BNS+ [16] Raef Bassily, Kobbi Nissim, Adam Smith, Thomas Steinke, Uri Stemmer, and Jonathan Ullman.  Algorithmic stability for adaptive data analysis.  In Proceedings of the forty-eighth annual ACM symposium on Theory of Computing, pages 1046–1059, 2016.
  * BST [14] Raef Bassily, Adam Smith, and Abhradeep Thakurta.  Private empirical risk minimization: Efficient algorithms and tight error bounds.  In 2014 IEEE 55th Annual Symposium on Foundations of Computer Science, pages 464–473. IEEE, 2014.
  * Chu [91] Alexander Mikhailovich Chudnov.  Game-theoretical problems of synthesis of signal generation and reception algorithms.  Problemy Peredachi Informatsii, 27(3):57–65, 1991.
  * CY [15] Yinzhi Cao and Junfeng Yang.  Towards making systems forget with machine unlearning.  In 2015 IEEE Symposium on Security and Privacy, pages 463–480. IEEE, 2015.
  * DMNS [06] Cynthia Dwork, Frank McSherry, Kobbi Nissim, and Adam Smith.  Calibrating noise to sensitivity in private data analysis.  In Theory of cryptography conference, pages 265–284. Springer, 2006\.
  * DMR [18] Luc Devroye, Abbas Mehrabian, and Tommy Reddad.  The total variation distance between high-dimensional gaussians.  arXiv preprint arXiv:1810.08693, 2018.
  * DR+ [14] Cynthia Dwork, Aaron Roth, et al.  The algorithmic foundations of differential privacy.  Foundations and Trends in Theoretical Computer Science, 9(3-4):211–407, 2014.
  * EKRR [19] Philip A Ernst, Wilfrid S Kendall, Gareth O Roberts, and Jeffrey S Rosenthal.  Mexit: Maximal un-coupling times for stochastic processes.  Stochastic Processes and their Applications, 129(2):355–380, 2019\.
  * GGHvdM [19] Chuan Guo, Tom Goldstein, Awni Hannun, and Laurens van der Maaten.  Certified data removal from machine learning models.  arXiv preprint arXiv:1911.03030, 2019.
  * GGVZ [19] Antonio Ginart, Melody Guan, Gregory Valiant, and James Y Zou.  Making ai forget you: Data deletion in machine learning.  In Advances in Neural Information Processing Systems, pages 3518–3531, 2019.
  * ISCZ [20] Zachary Izzo, Mary Anne Smart, Kamalika Chaudhuri, and James Zou.  Approximate data deletion from machine learning models: Algorithms and evaluations.  arXiv preprint arXiv:2002.10077, 2020.
  * Lan [12] Guanghui Lan.  An optimal method for stochastic composite optimization.  Mathematical Programming, 133(1-2):365–397, 2012.
  * LeC [98] Yann LeCun.  The mnist database of handwritten digits.  http://yann. lecun. com/exdb/mnist/, 1998.
  * LR+ [86] Torgny Lindvall, L Cris G Rogers, et al.  Coupling of multidimensional diffusions by reflection.  The Annals of Probability, 14(3):860–872, 1986.
  * Mir [17] Ilya Mironov.  Rényi differential privacy.  In 2017 IEEE 30th Computer Security Foundations Symposium (CSF), pages 263–275. IEEE, 2017.
  * NRSM [20] Seth Neel, Aaron Roth, and Saeed Sharifi-Malvajerdi.  Descent-to-delete: Gradient-based methods for machine unlearning.  arXiv preprint arXiv:2007.02923, 2020.
  * NY [83] Arkadij Semenovich Nemirovskij and David Borisovich Yudin.  Problem complexity and method efficiency in optimization.  Wiley-Interscience, 1983.
  * R+ [61] Alfréd Rényi et al.  On measures of entropy and information.  In Proceedings of the Fourth Berkeley Symposium on Mathematical Statistics and Probability, Volume 1: Contributions to the Theory of Statistics. The Regents of the University of California, 1961.
  * VEH [14] Tim Van Erven and Peter Harremos.  Rényi divergence and kullback-leibler divergence.  IEEE Transactions on Information Theory, 60(7):3797–3820, 2014\.
  * Ver [18] Roman Vershynin.  High-dimensional probability: An introduction with applications in data science, volume 47.  Cambridge university press, 2018.
  * Vil [08] Cédric Villani.  Optimal transport: old and new, volume 338.  Springer Science & Business Media, 2008.
  * Völ [16] Florian Völlering.  On maximal agreement couplings.  arXiv preprint arXiv:1608.01511, 2016.
  * Wal [77] Alastair J Walker.  An efficient method for generating discrete random variables with general distributions.  ACM Transactions on Mathematical Software (TOMS), 3(3):253–256, 1977.
  * Wik [21] Wikipedia.  Right to be forgotten — Wikipedia, the free encyclopedia.  http://en.wikipedia.org/w/index.php?title=Right%20to%20be%20forgotten&oldid=1007605238, 2021\.  [Online; accessed 23-February-2021].
  * WS [16] Blake E Woodworth and Nati Srebro.  Tight complexity bounds for optimizing composite objectives.  In D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 29. Curran Associates, Inc., 2016.

##  Appendix A Additional related work

We survey the works on machine unlearning - [11] 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. [5] 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. [16] 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
[18] 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^{2}d)$ 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 [23].

Our algorithm guarantees provable exact unlearning with probabilistic runtime
guarantees, whereas [23] 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
[23] 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 [23] 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, [23] (Theorem 3.4) achieve
an excess empirical risk of
$O\left({\frac{\sqrt{d}\sqrt{\operatorname{log}\left(1/\delta\right)}}{\epsilon
nk}}\right)^{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\left\\{{\frac{1}{\sqrt{\rho n}},\left({\frac{\sqrt{d}}{\rho
n}}\right)^{4/5}}\right\\}$ 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 kmT=k^{2}n\iff\rho=\frac{kn}{mT}$. Firstly, note that
as as long as $d\leq(\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 1, $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\geq(\rho n)^{3/4}$, our runtime is always smaller
than [23]: $kn$ as opposed to $k^{2}n$ gradient computations. Even with
$\rho=1$, our excess empirical risk is $\frac{1}{\sqrt{n}}$ and the excess
empirical risk of [23] is smaller than ours when
$\left({\frac{\sqrt{d}\sqrt{\operatorname{log}\left(1/\delta\right)}}{\epsilon
nk}}\right)^{2/5}\lesssim\frac{1}{\sqrt{n}}\iff\frac{\epsilon}{\sqrt{\operatorname{log}\left(1/\delta\right)}}\gtrsim\frac{\sqrt{d}n^{1/4}}{k}$.
In the second regime $d<(\rho n)^{3/4}$, we use Algorithm 5, wherein we have
$mT=\frac{(\rho n)^{2}}{d}$. This gives us $\rho=\frac{knd}{(\rho
n)^{2}}\iff\rho=\sqrt{\frac{kd}{n}}$, and our excess empirical risk is
$O\left({\left({\frac{\sqrt{d}}{\rho
n}}\right)^{4/5}}\right)=O\left({\frac{1}{(nk)^{2/5}}}\right)$. Therefore,
excess empirical risk of [23] is smaller than ours when
$\left({\frac{\sqrt{d}\sqrt{\operatorname{log}\left(1/\delta\right)}}{\epsilon
nk}}\right)^{2/5}\lesssim\frac{1}{(nk)^{2/5}}\iff\frac{\epsilon}{\sqrt{\operatorname{log}\left(1/\delta\right)}}\gtrsim
d$. We therefore have that unless $k$ is very large, the accuracy of [23] 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 [23], which will yield qualitatively
similar conclusions.

##  Appendix B Proofs for Section 5.1

###### Proof of Proposition 1.

We first show that Algorithm 1 is $\min\left\\{{1,\rho}\right\\}$-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
${\mathcal{A}}(S):=\widehat{\mathrm{w}}_{S}$ and
${\mathcal{A}}(S^{\prime}):=\widehat{\mathrm{w}}_{S^{\prime}}$ denote the
outputs of Algorithm 1 on $S$ and $S^{\prime}$ respectively. Since in
Algorithm 1, the randomness is only on indices, rather than actual data
points, say that $S=\left\\{{1,2,\ldots,n}\right\\}$. 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}=\left\\{{1,2,\ldots,n-1}\right\\}$. 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)=\begin{cases}\mu_{n-1,m}(b)&\text{ if }n\not\in b\\\
0&\text{otherwise}\end{cases}$. 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 $\mathrm{w}_{0}$ and other parameters,
Algorithm ${\mathcal{A}}(S)$ and ${\mathcal{A}}(S^{\prime})$ is the _same_
(deterministic) map from ${\mathbf{b}}=(b_{1},b_{2},\ldots,b_{T})$ where
$b_{j}\in[n]^{m}$ to ${\mathcal{W}}$. They only differ because of different
measures on the input space. Hence total variation distance between
${\mathcal{A}}(S)$ and ${\mathcal{A}}(S^{\prime})$ is just the total variation
distance between the push-forward measures
${\mathcal{A}}(S)_{\\#}\mu_{n,m}^{\otimes T}$ and
${\mathcal{A}}(S^{\prime})_{\\#}\mu_{n,m}^{{}^{\prime}\otimes T}$ which by
using the fact that ${\mathcal{A}}(S)\equiv{\mathcal{A}}(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,

$\displaystyle\text{TV}({\mathcal{A}}(S),{\mathcal{A}}(S^{\prime}))$
$\displaystyle\leq\text{TV}(\mu_{n,m}^{\otimes
T},\mu_{n,m}^{{}^{\prime}\otimes
T})=\sup_{{\mathbf{b}}\in\text{Pow}([n]^{m})^{T})}\left|{\mu_{n,m}^{\otimes
T}({\mathbf{b}})-\mu_{n,m}^{{}^{\prime}\otimes T}({\mathbf{b}})}\right|$
$\displaystyle=\mu_{n,m}^{\otimes T}\left({\mathbf{b}}\text{ such that at
least one }b_{j}\text{ contains }n\right)$ $\displaystyle\leq
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}\leq\frac{Tm}{n}$. Taking a uniform bound over all neighbouring
datasets $S^{\prime}$, we get that
$\sup_{\Delta(S,S^{\prime})=1}\text{TV}({\mathcal{A}}(S),{\mathcal{A}}(S^{\prime}))\leq\frac{Tm}{n}$.
By definition of TV distance, we trivially have that
$\sup_{\Delta(S,S^{\prime})=1}\text{TV}({\mathcal{A}}(S),{\mathcal{A}}(S^{\prime}))\leq
1$. Therefore, setting $m=\frac{\rho n}{T}$, we get the desired result that
the output of Algorithm 1 is $\min\left\\{{\rho,1}\right\\}$-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 $\mathrm{w}$, we have that

$\displaystyle\mathbb{E}_{b}\frac{\sum_{j\in b}\nabla
f(\mathrm{w},z_{j})}{m}=\sum_{\binom{n}{m}\text{ choices for
}b}\frac{\sum_{j\in b}\nabla
f(\mathrm{w},z_{j})}{m\binom{n}{m}}=\frac{\binom{n-1}{m-1}}{m\binom{n}{m}}\sum_{j=1}^{n}\nabla
f(\mathrm{w},z_{j})=\frac{\sum_{j=1}^{n}\nabla f(\mathrm{w},z_{j})}{n}$

where in the second equality, we use the observation that every
$\mathrm{z}_{j}$ appears in exactly $\binom{n-1}{m-1}$ terms over all choices
for $b$. We now bound its variance, denoted by a ${\mathcal{V}}^{2}$ by direct
computation.

$\displaystyle{\mathcal{V}}^{2}$
$\displaystyle=\mathbb{E}_{b}\left\|{\frac{\sum_{j\in b}\nabla
f(\mathrm{w},z_{j})}{m}-\mathbb{E}_{b}\left[\frac{\sum_{j\in b}\nabla
f(\mathrm{w},z_{j})}{m}\right]}\right\|^{2}$
$\displaystyle=\mathbb{E}_{b}\left\|{\frac{\sum_{j\in b}\nabla
f(\mathrm{w},z_{j})}{m}}\right\|^{2}-\left\|{\mathbb{E}_{b}\left[\frac{\sum_{j\in
b}\nabla f(\mathrm{w},z_{j})}{m}\right]}\right\|^{2}$
$\displaystyle=\sum_{{n\choose m}\text{ choices for }b}\frac{1}{{n\choose
m}}\frac{1}{m^{2}}\left\|{\sum_{j\in b}\nabla
f(\mathrm{w},z_{j})}\right\|^{2}-\left\|{\frac{\sum_{j=1}^{n}\nabla
f(\mathrm{w},z_{j})}{n}}\right\|^{2}$

In the first term, expanding the square and summing over all choices of $b$,
we get exactly ${n-1\choose m-1}$ terms of the form $\left\|{\nabla
f(\mathrm{w},\mathrm{z}_{j})}\right\|^{2}$ for $j=1$ to $n$, and ${n-2\choose
m-2}$ cross terms of the form $\left\langle{\nabla
f(\mathrm{w},\mathrm{z}_{i})},{\nabla
f(\mathrm{w},\mathrm{z}_{j})}\right\rangle$ for $i\neq 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

$\displaystyle{\mathcal{V}}^{2}$
$\displaystyle=\mathbb{E}_{b}\left\|{\frac{\sum_{j\in b}\nabla
f(\mathrm{w},z_{j})}{m}-\mathbb{E}_{b}\left[\frac{\sum_{j\in b}\nabla
f(\mathrm{w},z_{j})}{m}\right]}\right\|^{2}$
$\displaystyle=\left({\frac{{n-1\choose m-1}}{m^{2}{n\choose
m}}-\frac{1}{n^{2}}}\right)\sum_{j=1}^{n}\left\|{\nabla
f(\mathrm{w},\mathrm{z}_{j})}\right\|^{2}+\left({\frac{{n-2\choose
m-2}}{m^{2}{n\choose m}}-\frac{1}{n^{2}}}\right)\sum_{i,j=1,i\neq
j}^{n}\left\langle{\nabla f(\mathrm{w},\mathrm{z}_{i})},{\nabla
f(\mathrm{w},\mathrm{z}_{j})}\right\rangle$
$\displaystyle\leq\left({\frac{1}{mn}-\frac{1}{n^{2}}}\right)nG^{2}+\left|{\frac{m-1}{nm(n-1)}-\frac{1}{n^{2}}}\right|\sum_{i,j=1,i\neq
j}^{n}\left\|{\nabla f(\mathrm{w},\mathrm{z}_{i})}\right\|\left\|{\nabla
f(\mathrm{w},\mathrm{z}_{j})}\right\|$
$\displaystyle\leq\left({\frac{1}{m}-\frac{1}{n}}\right)G^{2}+\left|{\frac{m-1}{nm(n-1)}-\frac{1}{n^{2}}}\right|n(n-1)G^{2}$
$\displaystyle=\left({\frac{1}{m}-\frac{1}{n}}\right)G^{2}+\left|{\frac{m-1}{m}-\frac{(n-1)}{n}}\right|G^{2}$
$\displaystyle=\left({\frac{1}{m}-\frac{1}{n}}\right)G^{2}+\left|{\frac{1}{n}-\frac{1}{m}}\right|G^{2}$
$\displaystyle=2\left({\frac{1}{m}-\frac{1}{n}}\right)G^{2}\leq\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\leq 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 [3]) which when
using step size $\eta\leq\frac{1}{L}$ gives us

$\displaystyle\mathbb{E}{\widehat{F}_{S}(\widehat{\mathrm{w}}_{S})-\widehat{F}_{S}(\mathrm{w}^{*}_{S})}\leq
O\left({\frac{\eta{\mathcal{V}}^{2}}{(1-\eta L)}+\frac{D^{2}}{\eta T}}\right)$

Using step size $\eta\leq\frac{1}{2L}$, the right hand side simplifies to
$2\eta{\mathcal{V}}^{2}+\frac{D^{2}}{\eta
T}\leq\frac{4G^{2}\eta}{m}+\frac{D^{2}}{\eta T}=\frac{4G^{2}T\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\left\\{{\frac{1}{2L},\frac{D\sqrt{\rho n}}{GT}}\right\\}$ gives us

$\displaystyle\mathbb{E}{\widehat{F}_{S}(\widehat{\mathrm{w}}_{S})-\widehat{F}_{S}(\mathrm{w}^{*}_{S})}\leq
O\left({\frac{GD}{\sqrt{\rho n}}+\frac{D^{2}L}{T}}\right)$

Setting $T=\frac{DL\sqrt{\rho n}}{G}$ achieves the claimed result. ∎

###### Proof of Proposition 2.

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 [25]:

$\displaystyle D_{\alpha}(P\|Q)=\frac{1}{\alpha-1}\ln\left({\int
f_{P}(x)^{\alpha}f_{Q}(x)^{1-\alpha}dx}\right)$

Consider two neighbouring datasets $S=\left\\{{\mathrm{z}_{j}}\right\\}_{j}$
and $S^{\prime}=\left\\{{\mathrm{z}_{j}^{\prime}}\right\\}_{j}$ such that
$\Delta(S,S^{\prime})=1$, and let
$\left\\{{b_{t}^{\prime}}\right\\}_{t=1}^{T}$ and
$\left\\{{\mathrm{w}_{t}^{\prime}}\right\\}_{t=1}^{T}$ denote the mini-batch
indices and iterates of Algorithm 2 on dataset $S^{\prime}$ respectively. We
look at iteration $t$, and fix all the randomness before $t$ i.e. fix
$\mathrm{w}_{t}$ (and $\mathrm{w}_{t}^{\prime}$), as well as randomness in
sub-sampling mini-batch indices i.e. fix $b_{t}$. The $\alpha$-Rènyi
Divergence between $\mathrm{w}_{t+1}$ and $\mathrm{w}_{t+1}^{\prime}$ can be
bounded as,

$\displaystyle D_{\alpha}(\mathrm{w}_{t+1}\|\mathrm{w}_{t+1}^{\prime})$
$\displaystyle=D_{\alpha}\left({{\mathcal{P}}\left({\mathring{\mathrm{w}}_{t}-\eta\left({\frac{\sum_{j\in
b_{t}}\nabla
f(\mathring{\mathrm{w}}_{t},\mathrm{z}_{j})}{m}+\theta_{t}}\right)}\right)\Big{\|}{\mathcal{P}}\left({\mathring{\mathrm{w}}_{t}-\eta\left({\frac{\sum_{j\in
b_{t}}\nabla
f(\mathring{\mathrm{w}}_{t},\mathrm{z}_{j}^{\prime})}{m}+\theta_{t}}\right)}\right)}\right)$
$\displaystyle\leq
D_{\alpha}\left({\mathring{\mathrm{w}}_{t}-\eta\left({\frac{\sum_{j\in
b_{t}}\nabla
f(\mathring{\mathrm{w}}_{t},\mathrm{z}_{j})}{m}+\theta_{t}}\right)\Big{\|}\mathring{\mathrm{w}}_{t}-\eta\left({\frac{\sum_{j\in
b_{t}}\nabla
f(\mathring{\mathrm{w}}_{t},\mathrm{z}_{j}^{\prime})}{m}+\theta_{t}}\right)}\right)$
$\displaystyle\leq D_{\alpha}\left({\frac{\sum_{j\in b_{t}}\nabla
f(\mathring{\mathrm{w}}_{t},\mathrm{z}_{j})}{m}+\theta_{t}\Big{\|}\frac{\sum_{j\in
b_{t}}\nabla
f(\mathring{\mathrm{w}}_{t},\mathrm{z}_{j}^{\prime})}{m}+\theta_{t}}\right)\leq\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
$\left\|{\frac{\sum_{j\in b_{t}}\nabla
f(\mathring{\mathrm{w}}_{t},\mathrm{z}_{j})}{m}-\frac{\sum_{j\in b_{t}}\nabla
f(\mathring{\mathrm{w}}_{t},\mathrm{z}_{j}^{\prime})}{m}}\right\|^{2}\leq\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 [4], for $\alpha\leq 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 [4], when
considering $\alpha\leq 2$, the right hand side simplifies as
$\frac{1}{\alpha-1}\operatorname{log}\left(1+\frac{m^{2}}{n^{2}}\frac{\alpha(\alpha-1)}{2}4\left({\operatorname{exp}\left(\frac{8G^{2}}{m^{2}\sigma^{2}}\right)-1}\right)\right)\leq\frac{2\alpha
m^{2}}{n^{2}}\left({\operatorname{exp}\left(\frac{8G^{2}}{m^{2}\sigma^{2}}\right)-1}\right)\leq\frac{32\alpha
G^{2}}{n^{2}\sigma^{2}}$ where the last inequality use the numeric inequality
$\operatorname{exp}\left(x\right)\leq 1+2x$ when $x\leq 1.256$; this means
that we need the following condition $\frac{8G^{2}}{m^{2}\sigma^{2}}\leq
1.256$ \- we will revisit this condition later. The second point is that
Theorem $9$ in [4] holds integer $\alpha\geq 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 [26]). Therefore the result holds for all $\alpha\leq 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 [4] 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, [1]
for example. We now use adaptive sequential composition property of Rènyi
divergence (Proposition 1 in [22]) which linearly accumulates the divergence
across iterations, yielding that the Rènyi divergence between the iterates
$(\mathrm{w}_{1},\mathrm{w}_{2},\ldots,\mathrm{w}_{T})$ and
$(\mathrm{w}_{1}^{\prime},\mathrm{w}_{2}^{\prime},\ldots,\mathrm{w}_{T}^{\prime})$
is bounded as,
$D_{\alpha}((\mathrm{w}_{1},\mathrm{w}_{2},\ldots,\mathrm{w}_{T})\|(\mathrm{w}_{1}^{\prime},\mathrm{w}_{2}^{\prime},\ldots,\mathrm{w}_{T}^{\prime}))\leq\frac{64TG^{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 $\widehat{\mathrm{w}}_{S}$ and
$\widehat{\mathrm{w}}_{S}^{\prime}$. Moreover, this holds uniformly over all
neighbouring datasets $S^{\prime}$. We now use the result that
$\lim_{\alpha\rightarrow
1}D_{\alpha}(\widehat{\mathrm{w}}_{S}\|\widehat{\mathrm{w}}_{S}^{\prime})=D_{\text{KL}}(\widehat{\mathrm{w}}_{S}\|\widehat{\mathrm{w}}_{S}^{\prime})$
where $D_{\text{KL}}$ denotes the KL-divergence (see Theorem 5 in [26]). Hence
we get that
$D_{\text{KL}}(\widehat{\mathrm{w}}_{S}\|\widehat{\mathrm{w}}_{S}^{\prime})\leq\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}(\widehat{\mathrm{w}}_{S}\|\widehat{\mathrm{w}}_{S}^{\prime})\leq\sqrt{\frac{D_{\text{KL}}(\widehat{\mathrm{w}}_{S}\|\widehat{\mathrm{w}}_{S}^{\prime})}{2}}\leq\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}(\widehat{\mathrm{w}}_{S}\|\widehat{\mathrm{w}}_{S}^{\prime})\leq
1$; therefore setting $\sigma=\frac{8\sqrt{T}G}{n\rho}$ gives us that the
algorithm’s output is $\min\left\\{{\rho,1}\right\\}$ 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
Proposition 1 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,

$\displaystyle{\mathcal{V}}^{2}=\mathbb{E}{\left\|{\frac{\sum_{j\in
b_{t}}\nabla
f(\mathring{\mathrm{w}}_{t},\mathrm{z}_{j})}{m}+\theta_{t}-\nabla\widehat{F}_{S}(\mathring{\mathrm{w}}_{t})}\right\|}^{2}$
$\displaystyle=\mathbb{E}{\left\|{\frac{\sum_{j\in b_{t}}\nabla
f(\mathring{\mathrm{w}}_{t},\mathrm{z}_{j})}{m}-\nabla\widehat{F}_{S}(\mathring{\mathrm{w}}_{t})}\right\|}^{2}+\mathbb{E}{\left\|{\theta_{t}}\right\|}^{2}$
$\displaystyle\leq\frac{2G^{2}}{m}+\sigma^{2}d$

We now use Theorem 2 from [19] \- they use notation
$\left\\{{\beta_{t}}\right\\}_{t}$ and $\left\\{{\gamma_{t}}\right\\}_{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\leq\frac{1}{2L}$, and
appealing to Theorem 2 in [19], we get,

$\displaystyle\mathbb{E}{\widehat{F}(\widehat{\mathrm{w}}_{S})-\widehat{F}(\mathrm{w}^{*})}\leq
O\left({T\eta{\mathcal{V}}^{2}+\frac{D^{2}}{\eta T^{2}}}\right)=O\left({\eta
T\left({\frac{2G^{2}}{m}+\sigma^{2}d}\right)+\frac{D^{2}}{\eta T^{2}}}\right)$

Let $\tilde{G}^{2}=\frac{2G^{2}}{m}+\sigma^{2}d$, balancing the trade-off in
$\eta$ gives us $\eta=\frac{D}{\tilde{G}T^{3/2}}$. Therefore, setting
$\eta=\min\left\\{{\frac{1}{2L},\frac{D}{\tilde{G}T^{3/2}}}\right\\}$ gives us

$\displaystyle\mathbb{E}{\widehat{F}(\widehat{\mathrm{w}}_{S})-\widehat{F}(\mathrm{w}^{*})}$
$\displaystyle\leq
O\left({\frac{LD^{2}}{T^{2}}+\frac{\tilde{G}D}{\sqrt{T}}}\right)\leq
O\left({\frac{LD^{2}}{T^{2}}+\frac{GD}{\sqrt{Tm}}+\frac{\sigma\sqrt{d}D}{\sqrt{T}}}\right)$
$\displaystyle\leq
O\left({\frac{LD^{2}}{T^{2}}+\frac{GD}{\sqrt{Tm}}+\frac{GD\sqrt{d}}{n\rho}}\right)$

Finally, note that when using the amplification lemma, we arrived at the
condition $\frac{8G^{2}}{m^{2}\sigma^{2}}\leq 1.256$. Substituting
$\sigma=\frac{8\sqrt{T}G}{n\rho}$, this reduces to
$\frac{(n\rho)^{2}}{8m^{2}T}\leq 1.256\iff T\geq\frac{(n\rho)^{2}}{16m^{2}}$.

∎

###### Proof of Corollary 1.

We start with the result in Proposition 2, 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}}\geq\frac{LD^{2}}{T^{2}}\iff
m\leq\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\left({\frac{(\rho n)^{2}}{md},\sqrt{\frac{LD(\rho
n)}{G\sqrt{d}}}}\right)$ yields an expected excess empirical risk of
$O\left({\frac{GD\sqrt{d}}{n\rho}}\right)$.

We now look at the given condition $T\geq\frac{(n\rho)^{2}}{16m^{2}}$ given in
Proposition 2. We have set $T=\max\left({\frac{(\rho
n)^{2}}{md},\sqrt{\frac{LD(\rho n)}{G\sqrt{d}}}}\right)$, there we need to
ensure that $\frac{(\rho n)^{2}}{md}\geq\frac{(\rho n)^{2}}{16m^{2}}\iff
m\geq\frac{d}{16}$, as well as $\sqrt{\frac{LD(\rho
n)}{G\sqrt{d}}}\geq\frac{(\rho n)^{2}}{16m^{2}}\iff
m\geq\frac{1}{4}\left({\frac{(\rho n)^{3}G\sqrt{d}}{LD}}\right)^{1/4}$ \- this
recovers the condition
$m\geq\min\left\\{{\frac{d}{16},\frac{1}{4}\left({\frac{(\rho
n)^{3}G\sqrt{d}}{LD}}\right)^{1/4}}\right\\}$ in the Proposition statement.
Combining all the above arguments, we get that for any
$m\geq\min\left\\{{\frac{d}{16},\frac{1}{4}\left({\frac{(\rho
n)^{3}G\sqrt{d}}{LD}}\right)^{1/4}}\right\\}$ , setting
$T=\max\left\\{{\frac{(\rho n)^{2}}{md},\sqrt{\frac{LD(\rho
n)}{G\sqrt{d}}}}\right\\}$, yields an expected excess empirical risk of
$O\left({\frac{GD\sqrt{d}}{n\rho}}\right)$. ∎

######  Remark 4.

Note that in the above proof, if we use the stronger variance bound of
$2L^{2}\left({\frac{1}{m}-\frac{1}{n}}\right)$ from sub-sampling (derived in
the proof of Proposition 1), 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}}}$.

######  Proposition 5.

There exists neighbouring datasets $S$ and $S^{\prime}$ of $n$ points, and
smooth $G$-Lipshcitz convex functions $f$ and constraint set ${\mathcal{W}}$
such that the total variation distance between iterates produced by Algorithm
2 run on datasets $S$ and $S^{\prime}$, denoted by
$\left\\{{\mathrm{w}_{1},\mathrm{w}_{2},\ldots,\mathrm{w}_{T}}\right\\}$ and
$\left\\{{\mathrm{w}_{1}^{\prime},\mathrm{w}_{2}^{\prime},\ldots,\mathrm{w}_{T}^{\prime}}\right\\}$
respectively, is bounded as
$\text{TV}(\left({\mathrm{w}_{1},\mathrm{w}_{2},\ldots,\mathrm{w}_{T}}\right),\left({\mathrm{w}_{1}^{\prime},\mathrm{w}_{2}^{\prime},\ldots,\mathrm{w}_{T}^{\prime}}\right))\geq\min\left\\{{\Omega\left({\frac{G\sqrt{T}}{n\sigma}}\right),1}\right\\}$.

###### Proof of Proposition 5.

We first prove this without projection - let the constraint set
${\mathcal{W}}={\mathbb{R}}^{d}$, and so the projection ${\mathcal{P}}$ 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\mathrm{e}_{1}$ and that of $S^{\prime}$ be $G\mathrm{e}_{1}$, where
$\mathrm{e}_{1}$ is the first canonical basis vector. Let the function
$f(\mathrm{w},\mathrm{z})=\left\langle{\mathrm{w}},{\mathrm{z}}\right\rangle$.
The gradients are just data points $\mathrm{z}$, therefore gradients are $0$
on all but the differing points, wherein in the differing point in dataset
$S$, the gradient is a constant $-G\mathrm{e}_{1}$ and for dataset
$S^{\prime}$, it is $G\mathrm{e}_{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

$\displaystyle\text{TV}(\left({\mathrm{w}_{1},\mathrm{w}_{2},\ldots,\mathrm{w}_{T}}\right),\left({\mathrm{w}_{1}^{\prime},\mathrm{w}_{2}^{\prime},\ldots,\mathrm{w}_{T}^{\prime}}\right))\geq\text{TV}(\Psi\left({\mathrm{w}_{1},\mathrm{w}_{2},\ldots,\mathrm{w}_{T}}\right),\Psi\left({\mathrm{w}_{1}^{\prime},\mathrm{w}_{2}^{\prime},\ldots,\mathrm{w}_{T}^{\prime}}\right))=\text{TV}(\mathrm{w}_{T},\mathrm{w}_{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 $\mathrm{w}_{T}$ is
a mixture of multivariate Gaussians, all with variance
$T\sigma^{2}{\mathbb{I}}$ but with varying means:
$G\mathrm{e}_{1},2G\mathrm{e}_{1},\ldots,TG\mathrm{e}_{1}$, similarly for
$\mathrm{w}_{T}^{\prime}$. We denote the mixtures probabilities by $\pi_{i}$
where the $i^{th}$ conditional distribution, denoted by $\mathrm{w}_{T}^{i}$
and $\mathrm{w}_{T}^{\prime i}$ respectively, has means $iG\mathrm{e}_{1}$ and
$-iG\mathrm{e}_{1}$ respectively. Also, we denote the conditional probability
densities of the $i^{th}$ distribution by $\phi_{S}^{i}(\mathrm{w})$ and
$\phi_{S^{\prime}}^{i}(\mathrm{w})$ 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 $\left\\{{\mathrm{w}:\mathrm{w}_{1}\geq 0}\right\\}$. We can
therefore write the total variation distance as,

$\displaystyle\text{TV}(\mathrm{w}_{T}\|\mathrm{w}_{T}^{\prime})$
$\displaystyle=\frac{1}{2}\left\|{\phi_{S}(\mathrm{w})-\phi_{S^{\prime}}(\mathrm{w})}\right\|_{1}=\int_{\mathrm{w}_{1}\geq
0}\phi_{S}(\mathrm{w})-\phi_{S^{\prime}}(\mathrm{w})d\mathrm{w}=\int_{\mathrm{w}_{1}\geq
0}\sum_{i}\pi_{i}(\phi^{i}_{S}(\mathrm{w})-\phi^{j}_{S^{\prime}}(\mathrm{w}))d\mathrm{w}$
$\displaystyle=\sum_{i}\pi_{i}\int_{\mathrm{w}\geq
0}(\phi^{i}_{S}(\mathrm{w})-\phi^{j}_{S^{\prime}}(\mathrm{w}))d\mathrm{w}=\sum_{i}\pi_{i}\text{TV}(\mathrm{w}_{T}^{i},\mathrm{w}_{T}^{\prime
i})$
$\displaystyle\gtrsim\sum_{i}\pi_{i}\frac{2Gi}{m\sqrt{T}\sigma}=\frac{2G}{m\sqrt{T}\sigma}\mathbb{E}{i}=\frac{2G\sqrt{T}}{n\sigma}$

where in the inequality, we use the fact that $\mathrm{w}_{T}^{i}$ and
$\mathrm{w}_{T}^{\prime i}$ are Gaussians with means separated by $2G$, and
variance being $T\sigma^{2}{\mathbb{I}}$ and use the lower bound result on TV
between high-dimensional Gaussians [13]. 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
$\left\\{{\mathrm{w}:\mathrm{w}_{1}\geq 0}\right\\}$ 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\left\\{{\Omega\left({\frac{G\sqrt{T}}{n\rho}}\right),1}\right\\}$. ∎

##  Appendix C Proofs for Section 5.2

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 Del and Ins respectively. To elaborate, Del is a
(deterministic) map from $\left({[n]^{m}}\right)^{T}$ to
$\left({[n-1]^{m}}\right)^{T}$ and Ins is a map from
$\left({[n]^{m}}\right)^{T}$ to $\left({[n+1]^{m}}\right)^{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 3, the sub-sampled indices ${\mathbf{b}}$ are
drawn from a product distribution $\mu_{n,m}^{\otimes T}$ and in each
iteration of Algorithm 3 or Algorithm 4, the maps Del and 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.

###  C.1 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.

######  Claim 1 (Deletion).

For any set $b\in[n]^{m}$, we have that
$\mu^{\text{del}}_{n,m}(b)=\mu_{n-1,m}(b)$

###### Proof.

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 $\mathrm{z}$. 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
$\mathrm{z}\not\in b$. We therefore have,

$\displaystyle\mu_{n,m}^{\text{del}}(b)$
$\displaystyle=\mu_{n,m}(b|\left\\{{\mathrm{z}\not\in
b}\right\\})=\frac{\mu_{n,m}(b\cap\left\\{{\mathrm{z}\not\in
b}\right\\})}{\mu_{n,m}(\left\\{{\mathrm{z}\not\in b}\right\\})}$

By direct computation, $\mu_{n,m}(\left\\{{\mathrm{z}\not\in
b}\right\\})=1-\mu_{n,m}(\left\\{{\mathrm{z}\in
b}\right\\})=1-\frac{{n-1\choose m-1}}{{n\choose 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\left\\{{\mathrm{z}\not\in b}\right\\})=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\not\in b$. We have,

$\displaystyle\mu_{n,m}^{\text{del}}(b)$
$\displaystyle=\frac{\mu_{n,m}(b)}{\mu_{n,m}(\left\\{{\mathrm{z}\not\in
b}\right\\})}=\frac{1/{n\choose m}}{1-m/n}=\frac{n}{n-m}\frac{(n-m)!m!}{n!}$
$\displaystyle=\frac{(n-m-1)!m!}{(n-1)!}=\frac{1}{{n-1\choose
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.

######  Claim 2.

For any set $b\in[n+1]^{m}$, we have that
$\mu^{\text{ins}}_{n,m}(b)=\mu_{n+1,m}(b)$

###### Proof of 2.

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\left\\{{n+1}\right\\}\cup\left\\{{s_{i}}\right\\}|\left\\{{s_{i}}\right\\}\text{
replaced})\mu^{\text{ins}}_{n,m}(\left\\{{s_{i}}\right\\}\text{replaced}))=\mu_{n,m}(b\backslash\left\\{{n+1}\right\\}\cup\left\\{{s_{i}}\right\\}|\left\\{{s_{i}}\right\\})\nu(\left\\{{s_{i}}\right\\})=\frac{1}{{n\choose
m-1}}\frac{1}{(n+1-m)}$. We therefore have that

$\displaystyle\mu^{ins}_{n}(b)$
$\displaystyle=\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\choose
m-1}}\frac{1}{(n+1-m)}$ $\displaystyle=\frac{m}{n+1}\frac{1}{{n\choose
m-1}}=\frac{m(m-1)!(n-(m-1))!}{(n+1)n!}=\frac{1}{{n+1\choose
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

$\displaystyle\mu^{ins}_{n}(b)=\left({1-\frac{m}{n+1}}\right)\frac{1}{{n\choose
m}}=\frac{(n+1-m)(n-m)!m!}{n!(n+1)}=\frac{1}{{n+1\choose m}}=\mu_{n+1,m}(b)$

∎

##### Coupling.

We formally describe the coupling constructed by the unlearning Algorithm 3.
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
$\left\\{{c_{j}}\right\\}_{j=1}^{T}$, where
$c_{j}\sim\text{Bernoulli}\left({\frac{m}{n}}\right)$, 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 3 handles
insertions and deletions going from $S$ to $S^{\prime}$. We first show that
this is a valid coupling.

######  Claim 3.

For the coupling described above, for any $b$,

  1. 1.

$\mathbb{P}\left[{\mathbf{b}^{(1)}=b}\right]=\mu_{n}^{\otimes T}(b)$

  2. 2.

$\mathbb{P}\left[{\mathbf{b}^{(2)}=b}\right]=\mu_{n-1}^{\otimes T}(b)\ \
\text{(deletion)}$,
$\mathbb{P}\left[{\mathbf{b}^{(2)}=b}\right]=\mu_{n+1}^{\otimes T}(b)\ \ \
\text{(insertion)}$

###### Proof of 3.

Follows immediately from Claims 1 and 2 . ∎

We now show that the probability of disagreement under the above coupling is
upper bounded by $k$ times TV-stability parameter of Algorithm 1.

######  Claim 4.

For the $\rho$-TV stable Algorithm 1, under the coupling described above, the
following holds

$\displaystyle{\mathbb{P}}_{(\mathbf{b}^{(1)},\mathbf{b}^{(2)})}[\mathbf{b}^{(1)}\neq\mathbf{b}^{(2)}]\leq\rho$

###### Proof.

For deletion, we have,

$\displaystyle\mathbb{P}_{(\mathbf{b}^{(1)},\mathbf{b}^{(2)})}[\mathbf{b}^{(1)}\neq\mathbf{b}^{(2)}]=\mathbb{P}_{(\mathbf{b}^{(1)},\mathbf{b}^{(2)})}[\exists
j\in[T]:\mathbf{b}^{(1)}_{j}\neq\mathbf{b}^{(2)}_{j}]=\mathbb{P}_{\mathbf{b}}[\exists
j\in[T]:n\in\mathbf{b}_{j}]\leq\frac{Tm}{n}$

For insertion, we have

$\displaystyle\mathbb{P}_{(\mathbf{b}^{(1)},\mathbf{b}^{(2)})}[\mathbf{b}^{(1)}\neq\mathbf{b}^{(2)}]$
$\displaystyle={\mathbb{P}}_{{\mathbf{b}},\mathbf{c}}[\exists
j\in[T]:c_{j}=1]\leq\frac{Tm}{n}$

In Proposition 1, 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 Proposition 3.

###### Proof of Proposition 3.

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 Proposition
1, we embed the randomness for Algorithm 1 executed on $S$ and $S^{\prime}$
into a common probability space. Therefore, similar to the proof of
Proposition 1 given the datasets (and other parameters), Algorithm 1,
${\mathcal{A}}(S)$ is a deterministic map from sub-sampled indices
${\mathbf{b}}=(b_{1},b_{2},\ldots,\mathrm{b}_{T})$ to the model:
${\mathcal{A}}(S):{\mathbf{b}}\rightarrow{\mathcal{W}}$, 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 3. Hence it follows that the
output generated by Algorithm 1 has the same measure as
${\mathcal{A}}(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 4. Finally, from Remark 2, 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$. ∎

###  C.2 Unlearning for _noisy-m-A-SGD_

####  C.2.1 Coupling mini-batches

In this section, we show that Algorithm 4 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}$.

######  Claim 5.

For any set $b\in[n-1]^{m}$, we have that
$\mu^{\text{del}}_{n,m}(b)=\mu_{n-1,m}(b)$

###### Proof of 5.

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

$\displaystyle\mu^{\text{del}}_{n,m}(b)$
$\displaystyle=\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\left\\{{n}\right\\}$ 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\left\\{{n}\right\\})\nu(b_{i})=\frac{1}{{n-1\choose
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\choose m}}$.
Plugging these in, and with simple calculations, we have

$\displaystyle\mu^{\text{del}}_{n,m}(b)$
$\displaystyle=\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)$
$\displaystyle=\sum_{i=1}^{m}\mu_{n,m}(b\backslash
b_{i}\cup\left\\{{n}\right\\})\nu(b_{i})\frac{m}{n}+\frac{1}{{n-1\choose
m}}\left({1-\frac{m}{n}}\right)$ $\displaystyle=\frac{m}{{n-1\choose
m-1}}\frac{1}{n-1-(m-1)}\frac{m}{n}+\frac{1}{{n-1\choose
m}}\left({1-\frac{m}{n}}\right)$ $\displaystyle=\frac{1}{{n-1\choose
m}}+\frac{m}{n}\left({\frac{m}{{n-1\choose
m-1}(n-1-(m-1))}-\frac{1}{{n-1\choose m}}}\right)$
$\displaystyle=\frac{1}{{n-1\choose
m}}+\frac{m}{n}\left({\frac{m(m-1)!(n-1-(m-1))!}{(n-1)!(n-1-(m-1))}-\frac{1}{{n-1\choose
m}}}\right)$ $\displaystyle=\frac{1}{{n-1\choose
m}}+\frac{m}{n}\left({\frac{1}{{n-1\choose m}}-\frac{1}{{n-1\choose
m}}}\right)=\frac{1}{{n-1\choose m}}=\mu_{n-1,m}(b)\qed$

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}$.

######  Claim 6.