# Stochastic and Private Nonconvex Outlier-Robust PCA

Tyler Maunu  Chenyu Yu  Gilad Lerman

###### Abstract

We develop theoretically guaranteed stochastic methods for outlier-robust PCA.
Outlier-robust PCA seeks an underlying low-dimensional linear subspace from a
dataset that is corrupted with outliers. We are able to show that our methods,
which involve stochastic geodesic gradient descent over the Grassmannian
manifold, converge and recover an underlying subspace in various regimes
through the development of a novel convergence analysis. The main application
of this method is an effective differentially private algorithm for outlier-
robust PCA that uses a Gaussian noise mechanism within the stochastic gradient
method. Our results emphasize the advantages of the nonconvex methods over
another convex approach to solving this problem in the differentially private
setting. Experiments on synthetic and stylized data verify these results.

##  1 Introduction

_Outlier-robust PCA_ (ORPCA) involves the problem of robustly estimating an
underlying linear subspace from data in the presence of large amounts of
corrupted data. While many solutions have been proposed for this problem, some
particularly effective methods involve nonconvex energy minimization [26].
However, these methods require generic conditions on the full dataset and it
is not clear how they behave in the presence of _stochastic gradients_ , since
they typically require good initialization and control over where the iterates
lie.

This work develops a deeper understanding of how nonconvex methods for OR-PCA
interact with stochastic gradients. Past studies have mainly looked at
recovery limits of such methods [23], both in terms of percentages of
corrupted data as well as their associated statistics. In the current work, we
show that it is possible to extend the results to the stochastic setting while
maintaining robustness guarantees.

As an important application, we show that specific choices of stochastic
gradients lead to _differential privacy_. Private algorithms provide an
important way to gain insight from sensitive data. As a framework,
differential privacy has harkened in a new era in the study of privacy and its
interaction with data science and machine learning [14, 11]. To make an
algorithm differentially private, one typically incorporates some sort of
noise mechanism. This noise mechanism is applied to either the data itself or
within the algorithm to limit the influence any single point can have on the
output. In this paper, we focus on differentially private gradient descent
algorithms, which use noisy gradients at each iteration to achieve
differential privacy.

While differential privacy may be simple to include within an algorithm, it is
less straightforward to guarantee how accurate the algorithm will be.
Recently, there has been work on empirical risk minimization by differentially
private methods, which show that it is possible to achieve fast estimation and
optimization rates with differentially private algorithms [2, 30, 3]. Such
results typically focus on the convex setting, but some recent work has
studied such algorithms in the nonconvex setting as well [32]. While these
results are quite general, they do not capture the intricacies of the analysis
of robust methods. That is, especially in the setting of ORPCA, robust methods
are concerned with _recovery results_ , where under various conditions on a
corrupted dataset, an algorithm can still recover some underlying structure.
Especially in nonconvex recovery problems, it is not clear how the stochastic
nature of the private algorithms interacts with existing recovery guarantees.
Due to our generic guarantees for stochastic gradient methods, we are able to
guarantee recovery for a differentially private method.

###  1.1 Background

Suppose that we observe a dataset
$\mathcal{X}=\\{\bm{x}_{1},\dots,\bm{x}_{N}\\}\subset\mathbb{R}^{D}$. The
classical problem of principal component analysis (PCA) seeks the $r$
directions of maximum variance within this dataset, where $r$ is a parameter
chosen by the user. Equivalently, one can also try to find a linear subspace
that spans these directions. It is therefore convenient to encode PCA as a
problem over the Grassmannian manifold of $r$-dimensional linear subspaces in
$\mathbb{R}^{D}$, $G(D,r)$. Throughout the paper, we also consider the
optimization over orthogonal bases for $L\in G(D,r)$: each element of $G(D,r)$
can be spanned by the columns of a semiorthogonal matrix in
$O(D,r):=\\{\bm{V}\in\mathbb{R}^{D\times
r}:\bm{V}^{\top}\bm{V}=\bm{I}_{r}\\}$.

In this language, PCA solves the geometric problem

$\min_{\bm{V}\in
O(D,r)}\frac{1}{N}\sum_{i=1}^{N}\|(\bm{I}-\bm{V}\bm{V}^{\top})\bm{x}_{i}\|^{2},$
(PCA)

where $\bm{V}\bm{V}^{\top}$ is the orthogonal projection matrix onto
$\mathrm{Sp}(\bm{V})$. PCA thus finds the subspace which minimizes the sum of
squared distances between points and the subspace.

PCA is not outlier-robust due to the use of squared error. A typical way to
robustify it is to remove the square, which results in the following
formulation which we refer to as Grassmannian Least Absolute Deviations
(GLAD). :

$\displaystyle\min_{\bm{V}\in
O(D,r)}\frac{1}{N}\sum_{i=1}^{N}\|(\bm{I}-\bm{V}\bm{V}^{\top})\bm{x}_{i}\|=:F(\bm{V};\mathcal{X}).$
(GLAD)

Many methods have been proposed to solve this nonconvex and nonsmooth problem
and they are overviewed in Lerman and Maunu [23].

To give a high-level overview of our results, we will briefly discuss the two
areas that it straddles. First, the primary result of our analysis is
guarantees for a nonconvex, stochastic method for ORPCA. Typically, ORPCA
algorithms assume an inlier-outlier model,
$\mathcal{X}=\mathcal{X}_{\mathrm{in}}\cup\mathcal{X}_{\mathrm{out}}$, where
$\mathcal{X}_{\mathrm{in}}$ lie on a low-dimensional subspace $L_{\star}$, and
the outliers $\mathcal{X}_{\mathrm{out}}$ are corrupted to not lie on this
subspace. The goal is to recover $L_{\star}$, or $\bm{V}_{\star}\in O(D,r)$
such that $\mathrm{Sp}(\bm{V}_{\star})=L_{\star}$. For simplicity, we assume
that the data is centered, so that we search for a linear subspace. Throughout
the paper, we also make the simplifying assumption that $\mathcal{X}\subset
S^{D-1}$, where $S^{D-1}$ is the sphere in $\mathbb{R}^{D}$, so that the
function $F$ is 1-Lipschitz. This can be achieved by first normalizing all
points to the sphere, which has robustifying characteristics to adversarial
outliers [25].

Second, the important application of our results involves differential privacy
[13]. A randomized algorithm $\mathcal{A}$, which takes in an input $x$ and
gives back a random output, is $(\epsilon,\delta)$-differentially private if,
for all $\mathcal{S}\subseteq\text{Range}(\mathcal{A})$ and for all datasets
$x,y$ that only differ in at most one data point,
$P[\mathcal{A}(x)\in\mathcal{S}]\leq
e^{\epsilon}P[\mathcal{A}(y)\in\mathcal{S}]+\delta.$ Two common ways to make a
first-order algorithm private include adding noise to data or adding noise to
the gradients. In this work, since we study stochastic gradient methods for
outlier-robust PCA, the recovery guarantees we prove naturally extend to the
private setting.

###  1.2 Contributions

We derive the following results for ORPCA with large $N$:

  1. 1.

We present stochastic versions of the geodesic gradient descent (GGD)
algorithm, which results in the Noisy GGD (NGGD), Stochastic GGD (SGGD), and
Noisy Stochastic GGD (NSGGD) methods. We give theorems guaranteeing linear
convergence and subspace recovery by these three methods. Our results are the
first nonconvex convergence guarantees for stochastic gradient descent in the
least absolute deviations framework.

  2. 2.

With specifically chosen noise parameters, we demonstrate that these methods
are differentially private, and we refer to the resulting algorithms as dp-GGD
and dp-SGGD, respectively. We compare these private algorithms to convex
methods for differentially private outlier-robust PCA based on the REAPER
problem, (dp-REAP). In this setting, we extend past results on differentially
private convex empirical risk minimization to give subspace recovery
guarantees for the dp-REAP algorithms under generic conditions.

  3. 3.

By comparing our theoretical results for the differentially private methods,
we demonstrate a distinct advantage in the differentially private setting for
dp-(S)GGD over dp-REAP. The nonconvex dp-(S)GGD algorithm converges at a
linear rate while the convex dp-REAP methods converge at a sublinear rate,
meaning that one can obtain a much more accurate approximation to the
underlying subspace in less iterations. In terms of best approximations while
still maintaining privacy, we achieve approximation errors that are
$O(N^{-1})$ for the convex methods and errors on the order of
$O(2^{-N^{\tau}})$ for the nonconvex methods, where $\tau$ is some constant in
$(0,2)$ that depends on the statistics of the dataset.

  4. 4.

Experiments on synthetic and stylized data emphasize the theoretical results
of this paper. In particular, they demonstrate the advantage in terms of speed
and accuracy for the nonconvex methods, and in particular demonstrate distinct
advantages for the dp-SGGD method.

###  1.3 Review of Directly Related Work

For a comprehensive review of the many methods used for ORPCA, we direct the
reader to Lerman and Maunu [23]. Perhaps one of the most popular frameworks
for ORPCA uses least absolute deviations. Originating with the study of robust
orthogonal regression in Osborne and Watson [28], Späth and Watson [29], it
was considered for ORPCA in Ding et al. [10]. More recent studies by Zhang and
Lerman [35], Lerman et al. [24], Lerman and Maunu [22], Maunu et al. [26] have
demonstrated the considerable advantages of this program. This problem is
distinct from what is called Robust PCA (RPCA), which considers sparse
corruptions [6, 5].

The nonconvex method we propose is based on optimization on the Grassmannian
manifold [15]. Manifold optimization has recently been of great interest for
the machine learning community [34].

Differential privacy has become the preeminent way of protecting sensitive
data [13]. There has been a recent surge of work examining how differential
privacy affects the accuracy of various methods [2, 3]. Some recent work has
been devoted to considering differentially private methods for PCA [7, 20,
21].

###  1.4 Notation

We let $\sigma_{j}(\cdot)$ denote the $j$th singular value of a matrix. For
measuring subspace approximation, we use a distance metric on the
Grassmannian. A typical metric is
$d(L_{1},L_{2})=\sqrt{\sum_{j=1}^{r}\theta_{j}^{2}}$, where $\theta_{j}$ are
the principal angles between $L_{1}$ and $L_{2}$. For our later analysis of
the nonconvex method, for $\bm{V},\bm{V}^{\prime}\in O(D,r)$, which are bases
for two elements of $G(D,r)$, it is more convenient to work with the squared
metric
$d_{r}^{2}(\bm{V},\bm{V}^{\prime})=1-\sigma_{r}(\bm{V}^{\top}\bm{V}^{\prime})$,
which for subspaces that are close together is on the order of 1/2 times the
largest principal angle squared between $\mathrm{Sp}(\bm{V})$ and
$\mathrm{Sp}(\bm{V}^{\prime})$ (specifically, it is $1-\cos(\theta_{1})$). We
denote $B_{d_{r}^{2}}(\bm{V},\rho)$ to be the ball of radius $\rho$ with
respect to $d_{r}^{2}$.

##  2 Stochastic Algorithms to Minimize GLAD

In this paper, we propose to use stochastic gradient descent to directly
minimize (GLAD). This extends the existing framework for ORPCA studied by
Maunu et al. [26], where the authors proposed to use vanilla geodesic gradient
descent (GGD). Section 2.1 reviews the GGD method used to minimize (GLAD).
Then, Section 2.2 discusses modifications of this method to include noisy and
minibatch gradients, which result in stochastic GGD methods.

###  2.1 Geodesic Gradient Descent

One can directly optimize (GLAD) over the Grassmannian manifold using
geometric methods. Past algorithms that accomplish this with some theoretical
guarantees (despite the nonconvex setting) include IRLS [22] and GGD [26]. On
top of frequently being more accurate than their convex counterparts, these
methods are also faster than convex methods, since nonconvex methods work with
a $D\times r$ optimization variable rather than the typical $D\times D$
variable.

We briefly review GGD. Since $G(D,r)$ forms a Riemannian manifold, the
Riemannian gradient of the energy function in (GLAD) is

$\nabla
F(\bm{V};\mathcal{X})=\frac{1}{|\mathcal{X}|}\bm{Q}_{\bm{V}}\sum_{\bm{x}\in\mathcal{X}}\frac{\bm{x}\bm{x}^{\top}\bm{V}}{\|\bm{Q}_{\bm{V}}\bm{x}\|},$
(2.1)

where $\bm{V}\in O(D,r)$ is a matrix whose columns span $L$,
$\bm{Q}_{\bm{V}}=\bm{I}-\bm{V}\bm{V}^{\top}$ projects the gradient to the
tangent space of $G(D,r)$ and $|\mathcal{X}|$ denotes the number of points in
the set $\mathcal{X}$. Geodesic gradient descent (GGD) then takes the form
$\bm{V}_{k+1}=\operatorname{\mathsf{Exp}}_{\bm{V}_{k}}(-\eta_{k}\nabla
F(\bm{V}_{k};\mathcal{X}))$ (where $\operatorname{\mathsf{Exp}}$ is the
exponential map). For a complete discussion of this iteration and associated
concepts related to the geometry of $G(D,r)$ [15, 26].

###  2.2 Stochastic Geodesic Gradient Descent Methods

In terms of optimization, the main innovation in this work is to consider
stochastic gradient methods for (GLAD). One specific stochastic gradient one
may consider is the addition of Gaussian noise, which enhances privacy. To go
beyond this setting, we also consider stochasticity due to _minibatching_.
While the addition of stochastic gradients is a small modification of the
original GGD method, it is entirely nontrivial to extend convergence and
recovery analysis to the stochastic setting (see Section 3).

We first describe a version of GGD which uses noisy gradients. Let
$\bm{B}_{k}\in\mathbb{R}^{D\times r}$ whose entries are
i.i.d.$\sim\mathcal{N}(0,\sigma^{2})$. The noisy GGD (NGGD) iteration is given
by

${\bm{V}}_{k+1}=\mathcal{P}_{O(D,r)}({\bm{V}}_{k}-\eta_{k}(\nabla
F({\bm{V}}_{k};\mathcal{X})+\bm{B}_{k})),$ (NGGD)

where $\mathcal{P}_{O(D,r)}$ is the projection operator that solves
$\mathcal{P}_{O(D,r)}(\bm{A})=\mathop{\mathrm{argmin}}_{\bm{V}\in
O(D,r)}\|\bm{V}-\bm{A}\|^{2}$. This is an example of the orthogonal Procrustes
problem [19], and it can be solved via the SVD or polar decomposition [16].
This iteration is referred to as Noisy Geodesic Gradient Descent (NGGD).

We can also use stochastic estimates of $\nabla F({\bm{V}}_{k};\mathcal{X})$
to add further “noise” to the gradient. We call such a method noisy stochastic
geodesic gradient descent, NSGGD, which is defined by the iteration

$\tilde{\bm{V}}_{k+1}=\mathcal{P}_{O(D,r)}(\tilde{\bm{V}}_{k}-\eta_{k}(\bm{G}_{k}+\bm{B}_{k})).$
(NSGGD)

Here, $\bm{G}_{k}$ is an estimate of the gradient at $\tilde{\bm{V}}_{k}$.
When using minibatch stochastic gradients, we let
$\bm{G}_{k}=\mathsf{grad}F(\bm{V}_{k};\mathcal{B}_{k})$, where
$\mathcal{B}_{k}\subset\mathcal{X}$. We refer to the method with minibatch
stochastic gradients and zero noise as SGGD.

These methods have many potential applications. First, the minibatch SGGD
method allows for less per-iteration complexity than that of GGD, where SGGD
has complexity $O(|\mathcal{B}|Dd)$ per iteration and GGD has complexity of
$O(NDd)$ per iteration. Furthermore, the addition of noise allows for the
potential development of Langevin-like algorithms on the Grassmannian.
Finally, as we discuss later, when the noise has sufficiently large variance,
we can show that the resulting method is differentially private.

##  3 Theory

In the following sections we present our theoretical results for NGGD, SGGD,
and NSGGD. In particular, we prove convergence _and_ subspace recovery results
for these methods.

First, in Section 3.1, we recall a result from Maunu et al. [26], which shows
that PCA gives a good initial approximation to the underlying subspace with
high probability. After this, Section 3.2 gives an iteration complexity and
approximation result for NGGD. Then, Section 3.3 gives an iteration complexity
and approximation result for SGGD as well as a convergence and recovery
theorem for NSGGD. The proofs of convergence for these differentially private
methods require nontrivial extensions of the past proofs of convergence for
GGD seen in Maunu et al. [26]. After this, we finish in Section 3.4 by showing
how one can extend these approximation guarantees to achieve linear
convergence of the NGGD, SGGD, and NSGGD algorithms with a geometrically
diminishing step size scheme. For brevity, all proofs are left to the
Appendix.

The results in these sections represent the main theoretical innovation of
this work. Similar to the analysis of the deterministic GGD method, the
strategy is to prove, under a general condition called _stability_ , 1) good
initialization by some means, and 2) convergence of the nonconvex stochastic
gradient method.

###  3.1 Initialization by PCA

Our nonconvex methods require initialization in a sufficiently small
neighborhood of the true subspace spanned by $\bm{V}_{\star}$. To accomplish
this, we initialize NGGD, SGGD, and NSGGD using a PCA subspace. Later, in the
case of differentially private methods, we show that one can also initialize
with differentially private PCA. The main result for initialization follows.

######  Theorem 1 ([26]).

If

$\mathcal{S}^{\mathsf{PCA}}_{\gamma}(\mathcal{X}):=2\sin(\arccos(\gamma))\,\lambda_{r}(\bm{X}_{\mathrm{in}}\bm{X}_{\mathrm{in}}^{\top})-\|\bm{X}_{\mathrm{out}}\|_{2}^{2}>0,$
(3.1)

then $d_{r}^{2}(L_{PCA},L_{\star})<\gamma$.

###  3.2 Noisy GGD

Towards a complete theory for private, nonconvex robust subspace recovery, we
first prove an iteration complexity and approximation result for NGGD.
Following the analysis in Maunu et al. [26], the goal is to show that the
sequence $\sigma_{r}(\bm{V}_{\star}^{\top}\bm{V}_{k})$ forms a sequence that
rapidly increases with $k$.

For the convergence of GGD in Maunu et al. [26], the key idea is the
development of the stability statistic, which is defined as

$\displaystyle\mathcal{S}_{\gamma}(\mathcal{X})$
$\displaystyle=\gamma\lambda_{r}\left(\frac{1}{N}\sum_{\mathcal{X}_{\mathrm{in}}}\frac{\bm{x}\bm{x}^{\top}}{\|\bm{x}\|}\right)-\max_{\bm{V}\in
O(D,d)}\sigma_{1}\left(\nabla F(\bm{V};\mathcal{X}_{\mathrm{out}}))\right)$
$\displaystyle=\gamma\mathcal{P}(\mathcal{X}_{\mathrm{in}})-\mathcal{A}(\mathcal{X}_{\mathrm{out}}).$

Note that our parametrization of this statistic is slightly different from
that of Maunu et al. [26], where we use $\gamma$ instead of $\cos(\gamma)$.
Under _stability_ , or the assumption that
$\mathcal{S}_{\gamma}(\mathcal{X})>0$, Maunu et al. [26] prove local
convergence of GGD given initialization in $B(\bm{V}_{\star},\gamma)$. In the
following theorem, we prove convergence of NGGD when
$\mathcal{S}_{\gamma}(\mathcal{X})>0$ as long as $\bm{V}_{0}\in
B_{d_{r}^{2}}(\bm{V}_{\star},\gamma/2)$ – notice that the noisy method
requires some extra wiggle room.

######  Theorem 2.

Assume that $\mathcal{S}_{\gamma}(\mathcal{X})>0$, NGGD is initialized at
$\bm{V}_{0}\in B_{d_{r}^{2}}(\bm{V}_{\star},\gamma/2)$ with a constant step
size $\eta_{k}=s=c_{1}a/T^{\nu}$, $0.5<\nu<1$, and is run for $T$ iterations,
where

$\displaystyle
O(N^{2}\epsilon^{2})>T>\mathcal{F}_{1}(a/d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0}),\lambda),$

for $\mathcal{F}_{1}$ defined in (B.5). Then NGGD yields a final iterate
$\bm{V}_{T}\in{B(\bm{V}_{\star},a)}$ with probability at least $1-{2}\lambda$.

By Theorem 1, PCA initialization achieves the proper initialization with high
probability when the condition
$\mathcal{S}^{\mathsf{PCA}}_{\gamma}(\mathcal{X})>0$ holds. Theorem 2 states
that, effectively, as long as the number of iterations is
$>\mathcal{F}_{1}(c)$, the NGGD final iterate lies in
$B_{d_{r}^{2}}(\bm{V}_{\star},c\gamma)$ with high probability. We will show in
Section 3.4 how one can turn this into a linear convergence result.

###  3.3 Noisy Stochastic GGD

This section is mainly inspired by the analysis in Zhou et al. [36]. We first
present a novel analysis of minibatch SGGD for solving (GLAD).

To this end, we assume minibatches $\mathcal{X}^{k}$, $1\leq k\leq N_{B}$, of
size $B$ that are drawn from $\mathcal{X}$ with replacement. We can separate
each minibatch into inlier and outlier components
$\mathcal{X}_{\mathrm{in}}^{k}$ and $\mathcal{X}_{\mathrm{out}}^{k}$,
respectively. Much in the same way that one can analyze GGD and NGGD, we
analyze SGGD through the use of stability statistics. For SGGD, the main
difference is now each minibatch has an associated stability statistic,
$\mathcal{S}_{\gamma}(\mathcal{X}^{k})=\gamma\mathcal{P}(\mathcal{X}_{\mathrm{in}}^{k})-\mathcal{A}(\mathcal{X}_{\mathrm{out}}^{k})$.
As there are $N^{B}$ subsets $\mathcal{X}^{k}$, we get a range of stability
statistics, some of which are positive and some of which are negative. Now,
instead of assuming that $\mathcal{S}_{\gamma}(\mathcal{X})>0$, we assume that
for a minibatch selected uniformly at random from $\mathcal{X}$ with
replacement,

$\mathbb{E}\mathcal{S}_{\gamma}(\mathcal{X}^{k})=\mathcal{S}_{\gamma,\mathbb{E}}>0.$
(3.2)

######  Theorem 3.

Assume that $\mathcal{S}_{\gamma,\mathbb{E}}>0$, SGGD is initialized at
$\bm{V}_{0}\in B_{d_{r}^{2}}(\bm{V}_{\star},\gamma/2)$ with a constant step
size $\eta_{k}=s=c_{1}a/T^{\nu}$, $0.5<\nu<1$, and is run for $T$ iterations,
where

$O(N^{2}\epsilon^{2})>T>\mathcal{F}_{2}(a/d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0}),\lambda),$

for $\mathcal{F}_{2}$ defined in (B.10). Then SGGD yields a final iterate
$\bm{V}_{T}\in{B(\bm{V}_{\star},a)}$ with probability at least $1-{2}\lambda$.

To additionally prove convergence of NSGGD, we must also control the noise
throughout the iterations. This result essentially combines Theorems 2 and 3.
As before, this theorem states that in a number of iterations
$>\mathcal{F}_{3}(c)$, the NSGGD final iterate lies in
$B_{d_{r}^{2}}(\bm{V}_{\star},c\gamma)$ with high probability.

######  Theorem 4.

Assume that $\mathcal{S}_{\gamma,\mathbb{E}}>0$, NSGGD is initialized at
$\bm{V}_{0}\in B_{d_{r}^{2}}(\bm{V}_{\star},\gamma/2)$ with a constant step
size $\eta_{k}=s=c_{1}a/T^{\nu}$, $0.5<\nu<1$, and is run for $T$ iterations,
where

$O(N^{2}\epsilon^{2})>T>\mathcal{F}_{3}(a/d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0}),\lambda),$

for $\mathcal{F}_{3}$ defined in the Appendix. Then NSGGD yields a final
iterate $\bm{V}_{T}\in{B_{d_{r}^{2}}(\bm{V}_{\star},a)}$ with probability at
least $1-{4}\lambda$.

###  3.4 Linear Convergence Analysis

As we commented in the previous sections, in a constant number of iterations,
$\mathcal{F}_{1}(c)$ for NGGD, $\mathcal{F}_{2}(c)$ for SGGD, and
$\mathcal{F}_{3}(c)$ for NSGGD, the stochastic GGD algorithms converge to
$B_{d_{r}^{2}}(\bm{V}_{\star},c\gamma)$. Setting $c=1/2$, these methods have
yielded final estimates twice as close to $\bm{V}_{\star}$ in a constant
number of iterates.

Using this fact, the following theorem guarantees linear convergence of the
stochastic GGD algorithms. To accomplish this, we use a geometrically
diminishing step size. That is, we run the algorithm with a constant step size
$s$ for a sufficient number of iterations. Then, the algorithm is _restarted_
with a constant step size $cs$ for some fraction $c\in(0,1)$. This restarting
procedure is then repeated $R$ times. This is similar to the strategy used in
Maunu et al. [26] to prove linear convergence of GGD.

######  Theorem 5.

Suppose that one of the stochastic GGD algorithms is run for $R$ restarts and
$\mathcal{S}>0$, where $\mathcal{S}=\mathcal{S}_{\gamma}(\mathcal{X})$ for
NGGD and $\mathcal{S}=\mathcal{S}_{\gamma,\mathbb{E}}(\mathcal{X})$ for SGGD
and NSGGD. Suppose that in the first run of the algorithm (out of all the
restarts), the step size is $s=c_{1}aT_{1}^{-\nu}$, $0.5<\nu<1$ and the number
of iterations is
$T_{1}=O(\mathcal{F}_{j}(a/d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0})))$, where
$j=1$ for NGGD, $j=2$ for SGGD, and $j=3$ for $NSGGD$. Suppose further that
the step size for the $l$th restart is $s/2^{l-1}$ for
$T_{l}=\mathcal{F}_{j}(1/2)$ iterations, and
$\sum_{l=1}^{R}T_{l}=O(N^{2}\epsilon^{2})$. Then, with probability at least
$1-2R\lambda$ (or $1-4R\lambda$ for NSGGD), the output of the $R$th restart,
$\hat{\bm{V}}$, satisfies $\hat{\bm{V}}\in
B_{d_{r}^{2}}(\bm{V}_{\star},{a/2^{R}})$.

We see that this theorem guarantees an approximation that decreases at an
exponential rate over the number of restarts.

##  4 Application: Differential Privacy

In both the NGGD and NSGGD methods, if the noise variance is sufficiently
large, then the methods become differentially private. We guarantee the
privacy of these methods in the following theorem.

######  Theorem 6 (Differential Privacy of NGGD and NSGGD).

There exists a constant $c$ such that for any $\varepsilon<cT$, if
$\sigma^{2}\geq\frac{cT\log^{2}(1/\delta)}{\varepsilon^{2}N^{2}}$, then NGGD
run for $T$ iterations is $(\epsilon,\delta)$ differentially private. On the
other hand, if the batch size is $B$, there exist constants $c_{1}$ and
$c_{2}$ such that for any $\varepsilon<c_{1}q^{2}T$, if $\sigma^{2}\geq
c_{2}\frac{(B/N)^{2}T\log(1/\delta)}{\varepsilon^{2}N^{2}},$ then NSGGD run
for $T$ iterations is $(\epsilon,\delta)$ differentially private.

The proof of differential privacy for such stochastic first-order methods is
standard and follows Bassily et al. [2], Talwar et al. [30]. With the noise
variances as specified in Theorem 6, we refer to the NGGD algorithm as dp-GGD
and to NSGGD as dp-SGGD. When writing statements that apply to either dp-GGD
or dp-SGGD, we will refer to dp-(S)GGD.

In the following sections, we examine the implications of our recovery results
in the differentially private setting. First, Section 4.1 discusses how to
initialize NGGD and NSGGD in a private way. Then, Section 4.2 explains how the
results for NGGD and NSGGD translate to the differentially private setting.
Lastly, in Section 4.3 we present convex differentially private methods based
on REAPER [24], which gives an important baseline for subspace recovery based
on differentially private convex empirical risk minimization. In particular,
we extend convergence results for convex empirical risk minimization to the
case of the REAPER algorithm and show that these have implications for
subspace recovery.

###  4.1 PCA Initialization

Throughout the paper, we refer to dp-PCA as the output of the differentially
private PCA method of Jiang et al. [21]. Combining the previous result in
Theorem 1 with a result of Jiang et al. [21], we obtain the following theorem.

######  Theorem 7.

If $N$ is sufficiently large and
$\mathcal{S}^{\mathsf{PCA}}_{\gamma/2}(\mathcal{X})>0$, we have that the
output of dp-PCA, $\bm{V}_{dp-PCA}$, lies in
$B_{d_{r}^{2}}(\bm{V}_{\star},\gamma/2)$ with high probability.

###  4.2 Approximation for dp-(S)GGD

Notice that, in order for the conditions of Theorem 6 to be satisfied, we need
the total number of iterations to be bounded as $T=O(N^{2}\epsilon^{2})$. To
get a sense of the number of restarts we can take, we note that this implies
$\sum_{l=1}^{R}T_{l}=O(N^{2}\epsilon^{2})$. If we take $T_{l}=(\epsilon
N)^{\alpha}$ for $0<\alpha<2$ for all $l=1,\dots,R$, the conditions of Theorem
5 are satisfied once $N$ is sufficiently large. Therefore, we can take
$R=O((\epsilon N)^{2-\alpha})$, which yields a dp-(S)GGD estimator with
accuracy on the order of $1/2^{(\epsilon N)^{2-\alpha}}$, which decreases
exponentially in $N$. Taking this all together, we have the following
corollary of Theorems 5 and 6.

######  Corollary 1.

Running the dp-(S)GGD algorithm as in Theorem 5 for the maximum number of
restarts $R$ while still maintaining privacy in Theorem 6 yields a final
iterate $\hat{\bm{V}}$ such that
$d_{r}^{2}(\hat{\bm{V}},\bm{V}_{\star})=O(\gamma/2^{(\epsilon
N)^{2-\alpha}})$.

###  4.3 Differentially Private REAPER Algorithms

One could also attempt to relax (GLAD) and solve a surrogate convex problem
instead. A popular relaxation for this task is the REAPER relaxation of Lerman
et al. [24]. In this section, we present a simple differentially private
version of this method. We can directly apply existing empirical risk
minimization results to this problem [2, 3] to yield subspace recovery
guarantees. This will give us a baseline that demonstrates the superiority of
the nonconvex method.

The REAPER program [24] solves (GLAD) by relaxing the nonconvex constraints
that $\bm{P}_{L}$ is an orthoprojection:

$\displaystyle\min_{\bm{P}\in
H}\frac{1}{N}\sum_{i=1}^{N}\|(\bm{I}-\bm{P})\bm{x}_{i}\|,\
\mathcal{H}:=\\{\bm{P}:\bm{0}\preceq\bm{P}\preceq\bm{I},\
\operatorname{Tr}(\bm{P})=r\\}.$ (REAP)

This is a convex program, and so (REAP) can be solved by an array of standard
convex optimization algorithms. Since $\mathcal{X}\subset S^{D-1}$,
$G(\bm{P};\mathcal{X})=\frac{1}{N}\sum_{i=1}^{N}\|(\bm{I}-\bm{P})\bm{x}_{i}\|$
is 1-Lipschitz. Since the objective in (REAP) is not smooth, one must use
subgradient based methods [8]. We use the following subgradient of REAPER:

$\nabla
G(\bm{P};\mathcal{X})=-\sum_{\begin{subarray}{c}\bm{x}\in\mathcal{X}\\\
\|\bm{x}-\bm{P}\bm{x}\|>0\end{subarray}}\frac{(\bm{I}-\bm{P})\bm{x}\bm{x}^{T}+\bm{x}\bm{x}^{T}(\bm{I}-\bm{P})}{2\|\bm{x}-\bm{P}\bm{x}\|}.$
(4.1)

While Lerman et al. [24] proposes to solve this problem using an iteratively
reweighted least squares method, we instead opt to study first-order methods.
The first method we consider is gradient descent, and the second is a mirror
descent. To make these methods differentially private, we again use the
Gaussian mechanism and add noise to the gradient. Since past work has
demonstrated advantages for considering stochastic first-order methods when
making convex algorithms private [1, 3], we also give stochastic versions of
each algorithm. These convex optimization methods for the REAPER problem are
differentially private by the previous arguments of Bassily et al. [2], Talwar
et al. [30].

Since our primary focus is on the nonconvex method, and some nonprivate
versions of the convex methods were previously explored by Goes et al. [18],
we leave the exact formulation of these methods to the Appendix. In the
Appendix, we outline 4 differentially private algorithm for solving this
REAPER program: Differentially Private Gradient Descent (dp-GD-REAP),
Differentially Private Stochastic Gradient Descent (dp-SGD-REAP),
Differentially Private Mirror Descent (dp-MD-REAP), and Differentially Private
Stochastic Mirror Descent (dp-SMD-REAP).

Previous work on optimization with differential privacy has focused on
differentially private _empirical risk minimization_ [2, 30, 3]. In this
general set-up, one wishes to minimize the empirical surrogate for the
population loss. In the non-stochastic setting, we can use the main theorem of
Talwar et al. [30] for both dp-GD-REAP and dp-MD-REAP. Indeed, if one uses the
mirror map $\Psi(\cdot)=\|\cdot\|^{2}/2$, then the algorithm just becomes dp-
GD-REAP, whereas if one uses the negative von Neumann entropy, it yields dp-
MD-REAP. The following theorem gives our main approximation result for the
nonstochastic and stochastic dp-REAP algorithms. In both cases, we show that
the approximation error for these private methods is on $O(1/N)$, rather than
exponential like the dp-(S)GGD algorithms. In contrast to Bassily et al. [3],
this theorem does not resort to smoothing the cost function and instead uses
the optimization rate for subgradient descent.

######  Theorem 8.

Let $\mathsf{D}$ be the diameter of the constraint set $\mathcal{H}$. Then, if
dp-GD-REAP or dp-MD-REAP is run for $T=O(\epsilon^{2}N^{2})$ and yields the
estimator $\bar{\bm{P}}$, we have

$\mathbb{E}G(\bar{\bm{P}};\mathcal{X})-\min_{\bm{P}}G(\bm{P};\mathcal{X})\lesssim\frac{\mathsf{D}\log(N/\delta)}{\epsilon
N},$ (4.2)

where the expectation is taken over the randomness of the algorithm. On the
other hand, for dp-SGD-REAP and dp-SMD-REAP, if the noise variance is
$\sigma^{2}=c_{2}\frac{B^{2}T\log(1/\delta)}{\varepsilon^{2}N^{2}}$, then

$\displaystyle\mathbb{E}G(\bar{\bm{P}};\mathcal{X})-\min_{\bm{P}}G(\bm{P};\mathcal{X})$
$\displaystyle\leq\sqrt{1+c_{2}B^{2}\log(1/\delta)}\mathsf{D}\frac{1}{\epsilon
N},$

where the expectation is taken over the randomness of the algorithm and
$\mathsf{D}$ is the diameter of the constraint set with respect to the
appropriate geometry.

The proof for the nonstochastic methods is just Theorem 3.2 of Talwar et al.
[30], and the proof for the stochastic methods is given in Appendix B.4.

####  4.3.1 Implications for Subspace Recovery

We show that the approximate minimization guaranteed by Theorem 8 yields in a
generic setting approximate subspace recovery, or for REAPER, approximate
recovery of $\bm{P}_{\star}=\bm{P}_{L_{\star}}$.

We recall the following following _permeance_ , _alignment_ , and _stability_
statistics from Lerman et al. [24]:

$\displaystyle\mathcal{P}_{\mathsf{REAP}}$ $\displaystyle=\inf_{\bm{u}\in
L_{\star}\cap
S^{D-1}}\frac{1}{N}\sum_{\bm{x}\in\mathcal{X}_{\mathrm{in}}}|\bm{u}^{\top}\bm{P}_{\star}\bm{x}|,$
(4.3) $\displaystyle\mathcal{A}_{\mathsf{REAP}}$
$\displaystyle=\frac{1}{N}\|\bm{X}_{\mathrm{out}}\|\|\widetilde{\bm{Q}_{\star}\bm{X}_{\mathrm{out}}}\|,$
(4.4) $\displaystyle\mathcal{S}_{\mathsf{REAP}}$
$\displaystyle=\frac{\mathcal{P}}{4\sqrt{d}}-\mathcal{A}.$ (4.5)

Here, $\bm{Q}_{\star}=\bm{I}-\bm{P}_{\star}$ and the operator
$\widetilde{\cdot}$ normalizes the columns of
$\bm{Q}_{\star}\bm{X}_{\mathrm{out}}$ to the unit sphere. The permeance
measures how well spread the inliers are on the underlying subspace, the
alignment measures how aligned the outliers are orthogonal to $L_{\star}$, and
the stability is a tradeoff between these two terms. The result in Theorem 2.1
of Lerman et al. [24] states that if $\mathcal{S}_{\mathsf{REAP}}>0$, then
$\|\hat{\bm{P}}-\bm{P}_{\star}\|_{*}=0$, where $\|\cdot\|_{*}$ is the nuclear
or Schatten 1-norm. In other words, the REAPER program exactly recovers
$L_{\star}$ once $\mathcal{S}_{\mathsf{REAP}}>0$.

The following Theorem states the approximation result for the REAPER
algorithms of Section 4.3. In particular, it states that as $N$ increases and
the stability is bounded below, the distance between the REAPER subspace and
the true subspace goes to zero at a rate of $1/N$.

######  Theorem 9.

Suppose that $\mathcal{S}_{\mathsf{REAP}}>0$. If dp-MD-REAP or dp-GD-REAP are
run on the REAPER problem for $T=O(\epsilon^{2}N^{2})$ iterations, and if
$\hat{L}$ is the principal subspace of the output of one of these algorithms,
then

$\mathbb{E}d^{2}(\hat{L},L_{\star})\lesssim\frac{\pi}{2\mathcal{S}}\frac{\log(N/\delta)}{\epsilon
N}.$ (4.6)

On the other hand, if $\hat{L}$ is the principal subspace of the output of dp-
SMD-REAP or dp-SGD-REAP,

$\mathbb{E}d^{2}(\hat{L},L_{\star})\lesssim\frac{\pi}{2\mathcal{S}}\frac{\sqrt{1+c_{2}B^{2}\log^{2}(1/\delta)}}{\epsilon
N}.$ (4.7)

The GD, MD, SGD, and SMD algorithms differ in their respective constants.

###  4.4 Comparing Nonconvex and Convex Results

Notice that Theorem 9 and Corollary 1 use different distance metrics for
$G(D,r)$. It turns out that up to a factor of $r$, these are equivalent: for
an $r$-dimensional subspace, $d^{2}(L_{1},L_{2})\leq r\theta_{1}^{2}$, where
$\theta_{1}$ is the maximum principal angle between $L_{1}$ and $L_{2}$. On
the other hand,
$1-\sigma_{r}(\bm{V}_{1}^{\top}\bm{V}_{2})=1-\cos(\theta_{1})=O(\theta_{1}^{2})$
for small $\theta_{1}$. Up to a constant factor of $r$, the metrics
$d_{r}(\cdot,\cdot)$ and $d(\cdot,\cdot)$ are of the same order. Thus,
comparing the results of Theorem 9 to Corollary 1, we see that the nonconvex
methods have a distinct advantage in the private setting. That is, the convex
algorithm only achieves an approximation error of $O(N^{-1})$ while the
nonconvex methods achieve approximation errors that are $O(2^{-N^{\tau}})$

###  4.5 Stability and Privacy

We finish with a short discussion of the interaction between robustness and
privacy. Consider the stability result of the GGD algorithm, which states that
if $\mathcal{S}_{\gamma}(\mathcal{X})>0$, then GGD locally recovers the
underlying subspace $L_{\star}$. Notice that the robustness of the method
itself can yield privacy. This is stated in the following theorem.

######  Theorem 10.

Let $\mathcal{X}_{-i}$ be the dataset $\mathcal{X}$ with the $i$th datapoint
removed. Suppose that

$\displaystyle\mathcal{S}_{\gamma}(\mathcal{X}_{-i})>0,\
\mathcal{S}^{\mathsf{PCA}}_{\gamma}(\mathcal{X}_{-i})>0,\ i=1,\dots,N.$ (4.8)

Then, GGD with PCA initialization is differentially private.

While the condition of this theorem is hard to verify, it says that for
certain inlier-outlier datasets, one doesn’t even need to add noise to the GGD
algorithm, since _it is already private_. An in depth study of privacy is left
to future work, as the focus of this work is on the convergence of stochastic
GGD methods.

##  5 Differential Privacy Experiments

We performed experiments in order to demonstrate some of the predictions of
the substantial theory that was developed. The settings of our experiments
focus on differential privacy, but we emphasize that the results are more
general and similar experiments show the benefit of NGGD, SGGD, and NSGGD in
practice. Additional experiments are in the appendix.

###  5.1 Synthetic Experiments

We present two synthetic experiments in this section. The first tests the
convergence properties of the proposed algorithms for a setting with fixed
parameters. The second tests the methods over a range of sample sizes and
dimensions to look at their effect on subspace recovery. More comprehensive
experiments that demonstrate dependencies on other parameters are in the
supplemental material. All experiments were implemented on a PC with Intel
i7-9700 CPUs and 16GB RAM. Below, error refers to the distance between an
iterate and the underlying subspace,
$d^{2}(\mathrm{Sp}(\bm{V}_{k}),\mathrm{Sp}(\bm{V}_{\star}))$.

We tested the 6 proposed algorithms: dp-(S)GD-REAP, dp-(S)MD-REAP and
dp-(S)GGD. We set the step size for the 4 dp-REAP algorithms to be
$\eta_{k}=8/\sqrt{k}$. The step size for dp-GGD and dp-SGGD is
$\eta_{k}=1/2^{\lfloor k/50\rfloor}$. We use a fixed batch size
$B=\max\left(N\sqrt{\frac{\epsilon}{4T}},1\right)$ for both the convex and
nonconvex methods [3].

For both experiments, we randomly generate datasets according to the haystack
model, with Gaussian inliers and outliers, described in Lerman et al. [24].
Points are scaled to the sphere before running our methods.

In Figure 1, we plot the median and interquartile range of log-error versus
iteration for the six algorithms on 100 randomly generated sets. The fixed
model parameters are $r=2$, $D=20$, $N=2000$ and an inlier ratio 0.5. We set
the privacy parameters to be $\epsilon=0.8$ and $\delta=1/\sqrt{N}$. All
algorithms are run with $T=N$ iterations. We note that dp-(S)GGD converges
faster to the underlying subspace than dp-(S)GD-REAP and dp-MD-(S)REAP, since
its convergence rate is linear, unlike the sublinear rate for the convex
methods. Nevertheless, in the initial 600 iterations of dp-SGD-REAP and dp-
SMD-REAP, they converge at a faster rate than dp-SGGD (we also observe this
with the initial 100 iterations of the non-stochastic methods). If $D$ is not
large, it may be beneficial to initialize dp-(S)GGD with a corresponding dp-
REAP method instead of dp-PCA.

Figure 1: Convergence of the proposed algorithms with fixed parameters (see
main text). Each algorithm is repeated 100 times and the median log error with
a shaded region of interquartile range is plotted as a function of iterations
(left two figures) and time (right two figures).

Figure 2: Phase transition plot for the number of points, $N$, versus the
dimension, $D$. Each square represents the average log-distance between the
final iterate and the true subspace for each algorithm. The nonconvex methods
outperform the convex methods, and the stochastic method is able to perform
well for much larger $D$ due to the smaller noise required with stochastic
gradients.

For the second experiment, Figure 2 gives a phase transition plot of $N$ vs.
$D$. The data parameters are $r=2$, percentage of inliers is $0.5$, and the
total number of iterations of each algorithm is $T=2N$, and we set the privacy
parameters to be $\epsilon=0.8$ and $\delta=1/\sqrt{N}$. The step size for the
dp-(S)GD-REAP algorithms is $\eta_{k}=8/\sqrt{k}$, and the step size for
dp-(S)GGD is $\eta_{k}=1/2^{\lfloor k/50\rfloor}$. Each algorithm is run 50
times and we display the average of the log-errors of the final iterate. In
the non-stochastic case, the dp-GGD method outperforms the dp-GD-REAP method.
Furthermore, the stochastic versions take a smaller noise, and so the methods
are able to better approximate the underlying subspace for much larger $D$s.
Finally, as predicted by the theory, the approximations for dp-(S)GGD are much
more accurate than those for dp-(S)GD-REAP.

###  5.2 Stylized Application: Modified POPRES

To test on higher-dimensional data with some real characteristics, we create a
stylized dataset. It aims to imitate the Population Reference Sample (POPRES)
database extracted by Novembre et al. [27]. This highly private database
includes 3,192 European individuals with 500,568 alleles at SNP loci. Novembre
et al. [27] filtered SNPs and screened individuals to reduce the dataset to a
sample of $N=1,387$ individuals and $D=197,146$ SNPs. They applied PCA with
$r=2$ to the reduced data and demonstrated that the genetic information of the
selected sample correlates with a geographical map of Europe.

In view of our experience with the POPRES database, we find several issues
with directly using the procedure of Novembre et al. [27] when addressing the
machine learning community. First, POPRES is not publicly available. Second,
the suggested preprocessing of Novembre et al. [27] raises some questions
about the meaningful selection of reduced coordinates and individuals for
which a desired correlation with a given map can be demonstrated. Furthermore,
the reporting on the selection of individuals (supp. material of Novembre et
al. [27]) seems to reveal some private information.

In order to avoid these sensitive issues, we generated a stylized application
motivated by the work of Novembre et al. [27]. We used the publicly available
dataset provided on Github by the authors of Novembre et al. [27]. It was
obtained by applying (non-private) PCA with $r=20$ to their reduced data, so
the provided data matrix $\bm{X}$ has size $1387\times 20$. To simulate high-
dimensional SNP data and further privatize $\bm{X}$, we transform it as
follows: We chose $D=10,000$ and multiplied $\bm{X}$ by a random $20\times
10,000$ Gaussian orthogonal ensemble (GOE) matrix to get $\bm{X}^{\prime}$.
For outliers, we generated a $1,000\times 30$ random matrix of uniform i.i.d.
elements in $[-0.5,0.5]$ and multiplied this matrix by a random $30\times
10,000$ GOE matrix. We thresholded the inlier and outlier matrices to obtain
the three values -1, 0 and 1 to express alleles, which we then recode as 0, 1,
2 (see details in supplemental material). We concatenated the two matrices to
form an inlier-outlier $2,387\times 10,000$ matrix $\bm{Y}$ with elements in
$\\{0,1,2\\}$.

Figure 3: Recovered projections for the stylized POPRES dataset. Each
algorithm is run on a synthetically generated gene matrix which mimics the
original SNP data with $N=2387$ and $D=10000$. Out of these points, $1387$ lie
close to the underlying subspace, which recovers the shape of Europe, and 1000
outliers are generated to lie in the ambient space. We note that dp-SGGD is
able to recover the correct shape of Europe, unlike other algorithms.

Figure 3 demonstrates the application of dp-PCA [7], dp-SGGD, dp-GGD and dp-
SGD-REAP to $\bm{Y}$, and then plotting the projection of only the inliers
(which are also in $\bm{X}^{\prime}$). We note that dp-PCA, dp-GGD and dp-SGD-
REAP are unable to recover the target 2-dimensional subspace which indicates
the shape of Europe, whereas the embedding of dp-SGGD is relatively successful
in doing this. Additional figures and runtimes are included in the
supplemental material.

##  6 Conclusion

In this work, we initiate the first study of differentially private ORPCA
algorithms. Our results demonstrate the distinct advantages of taking a
nonconvex geometric approach to solving the ORPCA problem privately. In
particular, we show that the nonconvex dp-(S)GGD algorithm converges linearly
to the underlying subspace under a standard assumption of stability, in
contrast to the convex method that only converges sublinearly. The techniques
we use to guarantee the nonconvex dp-(S)GGD are interesting in their own right
because they are the first proofs of convergence for stochastic methods in
nonconvex formulation of least absolute deviations for ORPCA. Furthermore, our
experiments confirm our results and demonstrate the advantages of dp-GGD and
dp-SGGD. In fact, dp-SGGD seems to be superior to dp-GGD due to its ability to
use smaller noise in the Gaussian mechanism.

It would be interesting to further extend the NGGD results to studying the
mixing of Langevin dynamics.

There are directions to explore in future work. The following limitations are
of main interest to theoreticians. We only focus here on large $N$, where some
constants depend on $D$. It would be interesting to study the high-dimensional
regime; however, even current works on dp-PCA do not seem to apply to this
regime. There are also limitations due to the theoretical setting of ORPCA.
First, we consider here the common setting of inliers lying exactly on the
subspace and it would also be interesting to consider the interplay of noisy
inliers and differentially private subspace recovery. Second, we assume
centered data and search for linear subspaces and it will be good to justify
differentially private centering approaches or extend this work searching for
affine subspaces. Finally, robustness can enhance privacy [12], but we did not
explore this in the main text.

In addition, more practical limitations are as follows. First, we lack
experiments on real data, though we explained the difficulty of working with
and reporting results of the POPRES data. Second, while we theoretically
verify privacy, we do not yet have a good test to verify that the algorithms
are in fact private. Third, we do not know if our bounds are optimal, and it
would be good to tighten these results as well as prove lower bounds. Fourth,
the result for the nonconvex case is only local, and so it is unclear how the
methods perform in general settings. Fifth, the choice of parameters is not
sufficiently clear in the nonconvex case, and even in the convex case the
estimates are only approximate. At last, we require the strong assumption of
an inlier-outlier model, and it is not clear in general when data may meet
this assumption.

## References

  * Abadi et al. [2016] M. Abadi, H. B. McMahan, A. Chu, I. Mironov, L. Zhang, I. Goodfellow, and K. Talwar.  Deep learning with differential privacy.  In _Proceedings of the ACM Conference on Computer and Communications Security_ , 2016.  ISBN 9781450341394.  doi: 10.1145/2976749.2978318.
  * Bassily et al. [2014] R. Bassily, A. Smith, and A. 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.
  * Bassily et al. [2019] R. Bassily, V. Feldman, K. Talwar, and A. Guha Thakurta.  Private stochastic convex optimization with optimal rates.  _Advances in Neural Information Processing Systems_ , 32:11282–11291, 2019.
  * Bubeck et al. [2015] S. Bubeck et al.  Convex optimization: Algorithms and complexity.  _Foundations and Trends® in Machine Learning_ , 8(3-4):231–357, 2015.
  * Candès et al. [2011] E. J. Candès, X. Li, Y. Ma, and J. Wright.  Robust principal component analysis?  _Journal of the ACM (JACM)_ , 58(3):11, 2011.
  * Chandrasekaran et al. [2011] V. Chandrasekaran, S. Sanghavi, P. A. Parrilo, and A. S. Willsky.  Rank-sparsity incoherence for matrix decomposition.  _SIAM Journal on Optimization_ , 21(2):572–596, 2011.
  * Chaudhuri et al. [2013] K. Chaudhuri, A. D. Sarwate, and K. Sinha.  A near-optimal algorithm for differentially-private principal components.  _Journal of Machine Learning Research_ , 14, 2013.
  * Clarke [1990] F. H. Clarke.  _Optimization and nonsmooth analysis_.  SIAM, 1990.
  * Davis and Kahan [1970] C. Davis and W. M. Kahan.  The rotation of eigenvectors by a perturbation. iii.  _SIAM J. on Numerical Analysis_ , 7:1–46, 1970.
  * Ding et al. [2006] C. Ding, D. Zhou, X. He, and H. Zha.  $R_{1}$-PCA: rotational invariant $L_{1}$-norm principal component analysis for robust subspace factorization.  In _Proceedings of the 23rd International Conference on Machine Learning_ , pages 281–288. ACM, 2006.
  * Dwork [2008] C. Dwork.  Differential privacy: A survey of results.  In _International conference on theory and applications of models of computation_ , pages 1–19. Springer, 2008.
  * Dwork and Lei [2009] C. Dwork and J. Lei.  Differential privacy and robust statistics.  In _Proceedings of the forty-first annual ACM symposium on Theory of computing_ , pages 371–380, 2009.
  * Dwork and Roth [2013] C. Dwork and A. Roth.  The algorithmic foundations of differential privacy.  _Foundations and Trends in Theoretical Computer Science_ , 2013.  ISSN 15513068.  doi: 10.1561/0400000042.
  * Dwork et al. [2006] C. Dwork, F. McSherry, K. Nissim, and A. Smith.  Calibrating noise to sensitivity in private data analysis.  In _Theory of cryptography conference_ , pages 265–284. Springer, 2006.
  * Edelman et al. [1999] A. Edelman, T. A. Arias, and S. T. Smith.  The geometry of algorithms with orthogonality constraints.  _SIAM J. Matrix Anal. Appl._ , 20(2):303–353 (electronic), 1999.  ISSN 0895-4798.
  * Fan and Hoffman [1955] K. Fan and A. J. Hoffman.  Some metric inequalities in the space of matrices.  _Proceedings of the American Mathematical Society_ , 6(1):111–116, 1955.
  * Gawlik and Leok [2018] E. S. Gawlik and M. Leok.  High-order retractions on matrix manifolds using projected polynomials.  _SIAM Journal on Matrix Analysis and Applications_ , 39(2):801–828, 2018.
  * Goes et al. [2014] J. Goes, T. Zhang, R. Arora, and G. Lerman.  Robust stochastic principal component analysis.  In _Artificial Intelligence and Statistics_ , pages 266–274, 2014\.
  * Gower and Dijksterhuis [2004] J. C. Gower and G. B. Dijksterhuis.  _Procrustes problems_ , volume 30.  Oxford University Press on Demand, 2004.
  * Hardt and Price [2014] M. Hardt and E. Price.  The noisy power method: A meta algorithm with applications.  _Advances in neural information processing systems_ , 27:2861–2869, 2014.
  * Jiang et al. [2016] W. Jiang, C. Xie, and Z. Zhang.  Wishart mechanism for differentially private principal components analysis.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , page 1730–1736, 2016.
  * Lerman and Maunu [2018a] G. Lerman and T. Maunu.  Fast, robust and non-convex subspace recovery.  _Information and Inference: A Journal of the IMA_ , 7(2):277–336, 2018a.
  * Lerman and Maunu [2018b] G. Lerman and T. Maunu.  An overview of robust subspace recovery.  _Proceedings of the IEEE_ , 106(8):1380–1410, Aug 2018b.  ISSN 0018-9219.  doi: 10.1109/JPROC.2018.2853141.
  * Lerman et al. [2015] G. Lerman, M. B. McCoy, J. A. Tropp, and T. Zhang.  Robust computation of linear models by convex relaxation.  _Foundations of Computational Mathematics_ , 15(2):363–410, 2015.
  * Maunu and Lerman [2019] T. Maunu and G. Lerman.  Robust subspace recovery with adversarial outliers.  _arXiv preprint arXiv:1904.03275_ , 2019.
  * Maunu et al. [2019] T. Maunu, T. Zhang, and G. Lerman.  A well-tempered landscape for non-convex robust subspace recovery.  _JMLR_ , 2019.
  * Novembre et al. [2008] J. Novembre, T. Johnson, K. Bryc, Z. Kutalik, A. R. Boyko, A. Auton, A. Indap, K. S. King, S. Bergmann, M. R. Nelson, M. Stephens, and C. D. Bustamante.  Genes mirror geography within europe.  _Nature_ , 456(7218):98–101, Nov. 2008.  doi: 10.1038/nature07331.
  * Osborne and Watson [1985] M. R. Osborne and G. A. Watson.  An analysis of the total approximation problem in separable norms, and an algorithm for the total $l_{1}$ problem.  _SIAM Journal on Scientific and Statistical Computing_ , 6(2):410–424, 1985.
  * Späth and Watson [1987] H. Späth and G. A. Watson.  On orthogonal linear approximation.  _Numer. Math._ , 51:531–543, October 1987.
  * Talwar et al. [2014] K. Talwar, A. Thakurta, and L. Zhang.  Private empirical risk minimization beyond the worst case: The effect of the constraint set geometry.  _arXiv preprint arXiv:1411.5417_ , 2014.
  * Wang et al. [2017] D. Wang, M. Ye, and J. Xu.  Differentially private empirical risk minimization revisited: faster and more general.  In _Proceedings of the 31st International Conference on Neural Information Processing Systems_ , pages 2719–2728, 2017.
  * Wang et al. [2019] D. Wang, C. Chen, and J. Xu.  Differentially private empirical risk minimization with non-convex loss functions.  In _International Conference on Machine Learning_ , pages 6526–6535. PMLR, 2019.
  * Yu et al. [2015] Y. Yu, T. Wang, and R. J. Samworth.  A useful variant of the davis–kahan theorem for statisticians.  _Biometrika_ , 102(2):315–323, 2015.
  * Zhang and Sra [2016] H. Zhang and S. Sra.  First-order methods for geodesically convex optimization.  In _Conference on Learning Theory_ , pages 1617–1638. PMLR, 2016\.
  * Zhang and Lerman [2014] T. Zhang and G. Lerman.  A novel M-estimator for robust PCA.  _Journal of Machine Learning Research_ , 15(1):749–808, 2014.
  * Zhou et al. [2020] Z. Zhou, P. Mertikopoulos, N. Bambos, S. P. Boyd, and P. W. Glynn.  On the convergence of mirror descent beyond stochastic convex programming.  _SIAM Journal on Optimization_ , 30(1):687–716, 2020.

Supplementary Material

##  Appendix A REAPER Algorithms

In Algorithm 1, we give the dp-SGD-REAPER Algorithm. The full dp-GD-REAP
algorithm follows the same steps, but takes the full dataset as a batch at
each iteration, uses the noise variance $\sigma^{2}\geq
32T\log^{2}(T/\delta)/(\varepsilon^{2}N^{2})$ [30]. The projection step into
$\mathcal{H}$ that is in line 6 of Algorithm 1 can be implemented as the
water-filling procedure of Lerman et al. [24].

We can also use mirror descent to minimize the REAPER objective. This results
in the dp-MD-REAP algorithm, which we write in Algorithm 2. As in the previous
case, the full MD algorithm takes the full dataset as a batch at each
iteration and uses the noise variance $\sigma^{2}\geq
32T\log^{2}(T/\delta)/(\varepsilon^{2}N^{2})$ [30]. The mirror map is the von
Neumann entropy,
$\Psi(\cdot)=\frac{1}{4}(\operatorname{Tr}(\cdot\log(\cdot))+\log(D))$, and
this approach was used before in Goes et al. [18]. It turns out that the
Bregman projection for this choice of mirror map just corresponds to trace
renormalization.

Input: Dataset $\mathcal{X}$; subspace dimension $r$; step sizes $\eta_{k}$;
max number of step $T$; privacy parameters: $(\epsilon,\delta)$;
regularization parameter:$\lambda$; batch size $B$; noise variance
$\sigma^{2}=c_{2}\frac{(B/N)^{2}T\log(1/\delta)}{\varepsilon^{2}N^{2}}$.

1 $\bm{P}_{0}\leftarrow\bm{A}^{T}\bm{A}$, where $\bm{A}\in\mathbb{R}^{D\times
D}$ and $\bm{A}_{ij}\sim N(1,0.01)$;

2 for  _$k=1:T$_ do

3 Sample a batch $\mathcal{B}_{k}\subset\mathcal{X}$ of size $B$ uniformly
with replacement;

4 Sample $\bm{B}_{k}\in\mathbb{R}^{D\times D}$ such that
$\bm{B}_{k}=\bm{B}_{k}^{\top}$ and $\bm{B}_{k,ij}\sim N(0,\sigma^{2})$;

5 $\tilde{\nabla}F(\bm{P}_{k};\mathcal{B}_{k})=\nabla
F(\bm{P}_{k};\mathcal{B}_{k})+\bm{B}_{k}$ (using (4.1));

6 Update:
$\bm{P}_{k+1}:=\mathop{\mathrm{argmin}}_{\bm{P}\in\mathcal{H}}\Big{\|}\bm{P}-\Big{[}\left(\bm{P}_{k}-\eta_{k}\left(\tilde{\nabla}F(\bm{P}_{k};\mathcal{B}_{k})\right)\right)\Big{]}\Big{\|}$;

7

8 end for

Output: $\tilde{\bm{P}}\leftarrow\frac{1}{T}\sum_{k=1}^{T}\bm{P}_{k}$

Algorithm 1 dp-SGD-REAP

Input: Dataset $\mathcal{X}$; subspace dimension $r$; step sizes $\eta_{k}$;
max number of step $T$; privacy parameters: $(\epsilon,\delta)$;
regularization parameter:$\lambda$; batch size $B$; noise variance
$\sigma^{2}=c_{2}\frac{(B/N)^{2}T\log(1/\delta)}{\varepsilon^{2}N^{2}}$.

1 $\bm{P}_{0}\leftarrow\bm{A}^{T}\bm{A}$, where $\bm{A}\in\mathbb{R}^{D\times
D}$ and $\bm{A}_{ij}\sim N(1,0.01)$;

2 for  _$k=1:T$_ do

3 Sample a batch $\mathcal{B}_{k}\subset\mathcal{X}$ of size $B$ uniformly
with replacement;

4 Sample $\bm{B}_{k}\in\mathbb{R}^{D\times D}$ such that
$\bm{B}_{k}=\bm{B}_{k}^{\top}$ and $\bm{B}_{k,ij}\sim N(0,\sigma^{2})$;

5 $\tilde{\nabla}F(\bm{P}_{k};\mathcal{B}_{k})=\nabla
F(\bm{P}_{k};\mathcal{B}_{k})+\bm{B}_{k}$ (using (4.1));

6
$\bm{P}_{k+1}=\exp\left[\log(\bm{P}_{k})-\eta_{k}\tilde{\nabla}F(\bm{P}_{k};\mathcal{B}_{k})\right]$;

7 $\bm{P}_{k+1}=r\bm{P}_{k+1}/\operatorname{Tr}(\bm{P}_{k+1})$;

8

Output: $\tilde{\bm{P}}\leftarrow\frac{1}{T}\sum_{k=1}^{T}\bm{P}_{k}$

Algorithm 2 dp-SMD-REAP

##  Appendix B Supplemental Theory

###  B.1 Projection and Geodesic Gradient Descent

While the methods in NGGD and (NSGGD) are projected gradient methods, due to
Gawlik and Leok [17], these iteration very well approximate geodesics:

$\displaystyle\mathcal{P}_{O(D,r)}({\bm{V}}_{k}-$
$\displaystyle\eta_{k}(\nabla F({\bm{V}}_{k};\mathcal{X})+\bm{B}_{k}))=$ (B.1)
$\displaystyle\operatorname{\mathsf{Exp}}_{\bm{V}_{k}}(-\eta_{k}(\nabla
F(\bm{V}_{k};\mathcal{X})+\bm{B}_{k}))+O(\eta_{k}^{3}).$

###  B.2 Doob’s Maximial Inequality

For the convergence of the dp-GGD, SGGD, and dp-SGGD algorithms, we use the
following version of _Doob’s maximal inequality_.

######  Theorem 11 (Doob’s maximal inequality).

Suppose that $S_{t}=\sum_{i=1}^{t}X_{i}$ is a martingale with respect to the
filtered probability space $(\Omega,\mathcal{F},\\{\mathcal{F}_{t}\\},P)$,
then

$P\left(\max_{t=1,\dots,T}|S_{k}|\geq
x\right)\leq\frac{\mathbb{E}|S_{T}|^{2}}{x^{2}}.$ (B.2)

In the following proofs, the corresponding filtered probability space should
be apparent from context.

###  B.3 Differential Privacy of dp-GGD

###### Proof of Theorem 6.

The result for dp-GGD is just a rehash of the proof of Wang et al. [31,
Theorem 6.1] using the strong composition theorem. The result for dp-SGGD is
proven in Abadi et al. [1] using the moment accounting method. ∎

###  B.4 Convergence of Stochastic Mirror Descent

Here we prove Theorem 8.

###### Proof.

By a classic result that can be found, for example, in Bubeck et al. [4], if
the stochastic oracle is such that $\mathbb{E}\|(\nabla
G(\bm{P};\mathcal{B}))+\bm{B}\|^{2}\leq R^{2}$,

$\mathbb{E}G(\bar{\bm{P}})-\min_{\bm{P}}G(\bm{P})\leq
R\mathsf{D}\sqrt{\frac{2}{T}}.$

In this statement, the randomness is taken over the randomness of the
minibatches $\mathcal{B}$, as well as the randomness in the Gaussian noise
$\bm{B}$. We have

$\mathbb{E}\|(\nabla
G(\bm{P};\mathcal{B}))+\bm{B}\|^{2}\leq\mathbb{E}\|(\nabla
G(\bm{P};\mathcal{B}))\|^{2}+\mathbb{E}\|\bm{B}\|^{2}$

For the first, we have that $\mathbb{E}\|(\nabla
G(\bm{P};\mathcal{B}))\|^{2}\leq 1$. For the second,

$\displaystyle\mathbb{E}\|\bm{B}\|^{2}$
$\displaystyle=\sigma^{2}\mathbb{E}\|\bm{Z}\|^{2}\lesssim\sigma^{2}M$
$\displaystyle=c_{2}\frac{B^{2}T\log^{2}(1/\delta)}{\varepsilon^{2}N^{2}}M.$

Plugging in,

$\displaystyle\mathbb{E}G(\bar{\bm{P}})-\min_{\bm{P}}G(\bm{P})$
$\displaystyle\leq\sqrt{1+c_{2}\frac{B^{2}T\log^{2}(1/\delta)}{\varepsilon^{2}N^{2}}M}\mathsf{D}\sqrt{\frac{2}{T}}.$

If $T=O(N^{2}\epsilon^{2})$, then

$\displaystyle\mathbb{E}G(\bar{\bm{P}})-\min_{\bm{P}}G(\bm{P})$
$\displaystyle\leq\sqrt{1+c_{2}B^{2}\log^{2}(1/\delta)}\mathsf{D}\frac{1}{\epsilon
N}.$

∎

###  B.5 Proof of Theorem 9

Theorem 8 only guarantees an approximation to the minimizer of $F$ after $T$
iterations. To turn this then into a result of approximate recovery for REAPER
that we see in Theorem 9, we rely on the rate of ascent for the perturbed
objective by Lemma 2.3 of Lerman et al. [24].

######  Lemma 1 ([24], Lemma 2.3).

$G(\bm{P}_{\star}+\bm{\Delta})-G(\bm{P}_{\star})\geq\mathcal{S}\|\bm{\Delta}\|_{*}$

This states that the objective grows quickly when one is far from
$\bm{P}_{\star}$. Therefore, if the excess risk is on the order of $1/N$, then
the perturbation can also be bounded on the order of $1/N$. This is stated in
the following theorem.

######  Theorem 12.

Suppose that an algorithm to solve the REAPER problem yields a point
$\hat{\bm{P}}$ such that $F(\hat{\bm{P}})-G(\bm{P}_{\star})\leq\epsilon$. Let
$\hat{L}$ be the principal $r$-subspace of $\hat{\bm{P}}$. Then,

$d(\hat{L},L_{\star})\leq\frac{\pi}{2}\frac{\epsilon}{\mathcal{S}}$ (B.3)

###### Proof.

This follows from combining Lemma 1 with the Davis-Kahan $\sin\Theta$ Theorem
[9, 33]. ∎

The proof of Theorem 9 follows by combining Theorem 12 with Theorem 8.

###  B.6 dp-GGD Proofs

####  B.6.1 Proof of Theorem 7

Based on Jiang et al. [21], we have

$\displaystyle\cos(\theta_{1}(L_{dp-PCA},L_{\star}))$
$\displaystyle=d_{r}^{2}(\bm{V}_{dp-PCA},\bm{V}_{\star})\leq\|\bm{V}_{dp-
PCA}\bm{V}_{dp-PCA}^{\top}-\bm{V}_{\star}\bm{V}_{\star}^{\top}\|_{2}$
$\displaystyle\leq\|\bm{V}_{dp-PCA}\bm{V}_{dp-
PCA}^{\top}-\bm{V}_{PCA}\bm{V}_{PCA}^{\top}\|_{2}+\|\bm{V}_{PCA}\bm{V}_{PCA}^{\top}-\bm{V}_{\star}\bm{V}_{\star}^{\top}\|_{2}$
$\displaystyle\leq 2\sqrt{d}\|\bm{W}\|_{2}+\gamma/4$ $\displaystyle\leq
2\sqrt{d}O(d\log(d)/(N\epsilon))+\gamma/4,$

with high probability. For $N$ sufficiently large, we find
$d_{r}^{2}(\bm{V}_{dp-PCA},\bm{V}_{\star})<\gamma$ with high probability.

####  B.6.2 Proof of Theorem 2

###### Proof.

We can write

$\sigma_{r}(\bm{V}_{\star}^{\top}\bm{V}_{k+1})\geq\sigma_{r}(\bm{V}_{\star}^{\top}\bm{V}_{k})+\Big{(}\bm{\beta}_{1}^{\top}\bm{V}_{\star}^{\top}(-s(\mathsf{grad}F(\bm{V}_{k};\mathcal{X})+\bm{B}_{k})\bm{\beta}_{2}\Big{)}-c(\mathcal{X})s^{2}.$

Taking one minus both sides yields

$d_{r}^{2}(\bm{V}_{\star},\bm{V}_{k+1})\leq
d_{r}^{2}(\bm{V}_{\star}^{\top}\bm{V}_{k})+s\Big{(}\bm{\beta}_{1}^{\top}\bm{V}_{\star}^{\top}((\mathsf{grad}F(\bm{V}_{k};\mathcal{X})+\bm{B}_{k})\bm{\beta}_{2}\Big{)}+c(\mathcal{X})s^{2}.$

Using the fact that $\sin(\arccos(x))=\sqrt{1-x^{2}}\geq\sqrt{1-x}$, $x\geq
0$, stability implies that

$s\Big{(}\bm{\beta}_{1}^{\top}\bm{V}_{\star}^{\top}((\mathsf{grad}F(\bm{V}_{k};\mathcal{X})+\bm{B}_{k})\bm{\beta}_{2}\Big{)}\leq-s\mathcal{S}_{\gamma}(\mathcal{X})d_{r}(\bm{V}_{\star},\bm{V}_{k})+s\bm{\beta}_{1}^{\top}\bm{V}_{\star}^{\top}\bm{B}_{k}\bm{\beta}_{2}$

Thus

$d_{r}^{2}(\bm{V}_{\star},\bm{V}_{k+1})\leq\Big{(}1-\frac{s\mathcal{S}_{\gamma}(\mathcal{X})}{d_{r}(\bm{V}_{\star}^{\top}\bm{V}_{k})}\Big{)}d_{r}^{2}(\bm{V}_{\star},\bm{V}_{k})+s\bm{\beta}_{1}^{\top}\bm{V}_{\star}^{\top}\bm{B}_{k}\bm{\beta}_{2}+c(\mathcal{X})s^{2}.$
(B.4)

Let $m_{T}=\max_{j=1,\dots,T}d_{r}(\bm{V}_{\star},\bm{V}_{j})$. We can
iteratively apply this to yield

$\displaystyle
d_{r}^{2}(\bm{V}_{\star},\bm{V}_{T})\leq(1-s\mathcal{S}_{\gamma}(\mathcal{X})/m_{T})^{T}d_{r}^{2}(\bm{V}_{\star},\bm{V}_{T})+$
$\displaystyle\sum_{j=1}^{T}s(1-s\mathcal{S}_{\gamma}(\mathcal{X})/m_{T})^{T-j}\bm{\beta}_{1}^{\top}\bm{V}_{\star}^{\top}\bm{B}_{k}\bm{\beta}_{2}+$
$\displaystyle\sum_{j=1}^{T}(1-s\mathcal{S}_{\gamma}(\mathcal{X})/m_{T})^{T-j}c(\mathcal{X})s^{2}.$

Bounding the maximum: We proceed by first bounding $m_{T}$ in terms of
$d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0})$ and other quantities. To do this, first
notice that this amounts to bounding $d_{r}^{2}(\bm{V}_{\star},\bm{V}_{k})$
for $k=1,\dots,T$. We can use (B.4) along with $\sin(\arccos(x))\geq 1-x$ to
write

$d_{r}^{2}(\bm{V}_{\star},\bm{V}_{k+1})\leq\Big{(}1-s\mathcal{S}_{\gamma}(\mathcal{X})\Big{)}d_{r}^{2}(\bm{V}_{\star},\bm{V}_{k})+s\bm{\beta}_{1}^{\top}\bm{V}_{\star}^{\top}\bm{B}_{k}\bm{\beta}_{2}+c(\mathcal{X})s^{2}.$

We proceed by applying Doob’s maximal inequality

$\Pr\Big{(}\max_{k=1,\dots,T}\Big{|}\sum_{j=1}^{k}s(1-s\mathcal{S}_{\gamma}(\mathcal{X}))^{k-j}\bm{\beta}_{1}^{j\top}\bm{V}_{\star}^{\top}\bm{B}_{j}\bm{\beta}_{2}^{j}\Big{|}>\epsilon\Big{)}\\\
\leq\Pr\Big{(}\max_{k=1,\dots,T}\Big{|}\sum_{j=1}^{k}s\bm{\beta}_{1}^{j\top}\bm{V}_{\star}^{\top}\bm{B}_{j}\bm{\beta}_{2}^{j}\Big{|}>\epsilon\Big{)}\leq\frac{s^{2}T\sigma^{2}(\sqrt{D}+\sqrt{d})^{2}}{\epsilon^{2}}.$

Choosing $\epsilon=\frac{s\sqrt{T}\sigma(\sqrt{D}+\sqrt{d})}{\sqrt{\lambda}}$
yields

$\displaystyle\Pr\Big{(}\max_{k=1,\dots,T}\Big{|}\sum_{j=1}^{k}s(1-s\mathcal{S}_{\gamma}(\mathcal{X}))^{k-j}\bm{\beta}_{1}^{j\top}\bm{V}_{\star}^{\top}\bm{B}_{j}\bm{\beta}_{2}^{j}\Big{|}>\epsilon\Big{)}$
$\displaystyle\leq\lambda.$

With probability at least $1-\lambda$, we thus have that for all
$k=1,\dots,T$,

$d_{r}^{2}(\bm{V}_{\star},\bm{V}_{k})\leq\Big{(}1-s\mathcal{S}_{\gamma}(\mathcal{X})\Big{)}^{k}d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0})+\frac{s\sqrt{k}\sigma(\sqrt{D}+\sqrt{d})}{\sqrt{\lambda}}+\frac{c(\mathcal{X})sk}{\mathcal{S}_{\gamma}(\mathcal{X})}.$

Thus, if

$\frac{s\sqrt{T}\sigma(\sqrt{D}+\sqrt{d})}{\sqrt{\lambda}}<\frac{d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0})}{2}$

and

$\frac{s^{2}Tc(\mathcal{X})}{\mathcal{S}_{\gamma}(\mathcal{X})}<\frac{d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0})}{2},$

then $m_{T}<2d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0})$. In particular, for
$s=c_{1}aT^{-\nu}$, these are satisfied if

$T>\max\Big{(}\Big{[}\frac{2c_{1}a^{2}c(\mathcal{X})}{\mathcal{S}_{\gamma}(\mathcal{X})d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0})}\Big{]}^{1/(2\nu-1)},\Big{[}\frac{2a\sigma(\sqrt{D}+\sqrt{d})}{\sqrt{\lambda}d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0})}\Big{]}^{2/(2\nu-1)}\Big{)}.$

Since $a>a^{2}$, a sufficient condition is

$T>\max\Big{(}\Big{[}\frac{2c_{1}ac(\mathcal{X})}{\mathcal{S}_{\gamma}(\mathcal{X})d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0})}\Big{]}^{1/(2\nu-1)},\Big{[}\frac{2a\sigma(\sqrt{D}+\sqrt{d})}{\sqrt{\lambda}d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0})}\Big{]}^{2/(2\nu-1)}\Big{)}=:\mathcal{F}_{1}(a/d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0}),\lambda).$
(B.5)

Notice that this is a function of $a/d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0})$ and
$\lambda$.

Choosing
$s=\frac{Ca}{T^{\nu}}\leq\frac{Cd_{r}^{2}(\bm{V}_{\star},\bm{V}_{0})}{2\sqrt{T}}$
then yields that $m_{T}<2d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0})$. In particular,
as long as we initialize in $B_{d_{r}^{2}}(\bm{V}_{\star},\gamma/2)$, we see
that stability holds throughout all iterations with probability at least
$1-\lambda$.

Bounding the $T$th iterate: Let $m_{0}=2d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0})$.
From here the proof is straightforward: the first term geometrically
decreases. The second can be bounded with Doob’s maximal inequality with high
probability and uses the fact that

$\sum_{j}(1-s\mathcal{S}_{\gamma}(\mathcal{X})/m_{0})^{2(T-j)}\leq\frac{1}{1-(1-s\mathcal{S}_{\gamma}(\mathcal{X})/m_{0})^{2}}=\frac{m_{0}^{2}}{2s\mathcal{S}_{\gamma}(\mathcal{X})m_{0}-(s\mathcal{S}_{\gamma}(\mathcal{X}))^{2}},$

which is independent of $T$. More specifically, Doob’s maximal inequality
yields

$\Pr\Big{(}\inf_{1\leq k\leq
T}\Big{|}\sum_{j=1}^{k}s(1-s\mathcal{S}_{\gamma}(\mathcal{X})/m_{0})^{T-j}\bm{\beta}_{1}^{j\top}\bm{V}_{\star}^{\top}\bm{B}_{j}\bm{\beta}_{2}^{j}\Big{|}>\epsilon\Big{)}\leq\\\
\frac{\mathbb{E}\Big{(}\sum_{j=1}^{T}s^{2}(1-s\mathcal{S}_{\gamma}(\mathcal{X})/m_{0})^{2(T-j)}\Big{[}\bm{\beta}_{1}^{j\top}\bm{V}_{\star}^{\top}\bm{B}_{j}\bm{\beta}_{2}^{j}\Big{]}^{2}\Big{)}}{\epsilon^{2}}.$

We can upper bound

$\displaystyle\mathbb{E}(\bm{\beta}_{1}^{j\top}\bm{V}_{\star}^{\top}\bm{B}_{j}\bm{\beta}_{2}^{j})^{2}$
$\displaystyle=\mathbb{E}(\bm{\alpha}_{j}^{\top}\bm{B}_{j}\bm{\beta}_{2}^{j})^{2}\leq\sigma^{2}\mathbb{E}\sigma_{1}(\bm{Z}_{j})^{2}$
$\displaystyle\lesssim\sigma^{2}(\sqrt{D}+\sqrt{d})^{2}$

In any case, this implies that

$\Pr\Big{(}\inf_{1\leq k\leq
T}\Big{|}\sum_{j=1}^{k}s(1-s\mathcal{S}_{\gamma}(\mathcal{X})/m_{0})^{T-j}\bm{\beta}_{1}^{j\top}\bm{V}_{\star}^{\top}\bm{B}_{j}\bm{\beta}_{2}^{j}\Big{|}\leq\\\
C\frac{\sigma(\sqrt{D}+\sqrt{d})\sqrt{\sum_{j=1}^{T}s^{2}(1-s\mathcal{S}_{\gamma}(\mathcal{X})/m_{0})^{2(T-j)}}}{\sqrt{\lambda}}\Big{)}\leq\lambda,$

or

$\Pr\Big{(}\inf_{1\leq k\leq
T}\Big{|}\sum_{j=1}^{k}s(1-s\mathcal{S}_{\gamma}(\mathcal{X})/m_{0})^{T-j}\bm{\beta}_{1}^{j\top}\bm{V}_{\star}^{\top}\bm{B}_{j}\bm{\beta}_{2}^{j}\Big{|}\leq\\\
Cs\frac{\sigma(\sqrt{D}+\sqrt{d})}{\sqrt{\lambda}}\sqrt{\frac{m_{0}^{2}}{2s\mathcal{S}_{\gamma}(\mathcal{X})m_{0}-(s\mathcal{S}_{\gamma}(\mathcal{X}))^{2}}}\Big{)}\leq\lambda.$

The last term uses the fact that

$\sum_{j=1}^{T}(1-s\mathcal{S}_{\gamma}(\mathcal{X})/m_{0})^{T-j}\leq\frac{m_{0}}{s\mathcal{S}_{\gamma}(\mathcal{X})}$

Putting these together, we find with probability at least $1-2\lambda$,

$d_{r}^{2}(\bm{V}_{\star},\bm{V}_{T})\leq(1-s\mathcal{S}_{\gamma}(\mathcal{X})/m_{0})^{T}d_{r}^{2}(\bm{V}_{\star}^{\top}\bm{V}_{0})\\\
+s\Big{[}C\frac{\sigma(\sqrt{D}+\sqrt{d})}{\sqrt{\lambda}}\frac{m_{0}}{\sqrt{2s\mathcal{S}_{\gamma}(\mathcal{X})m_{0}-(s\mathcal{S}_{\gamma}(\mathcal{X}))^{2}}}+\frac{m_{0}c(\mathcal{X})}{\mathcal{S}_{\gamma}(\mathcal{X})}\Big{]}.$

Now, if $T$ is sufficiently large so that $s=\frac{Ca}{T^{\nu}}$ satisfies

$\displaystyle s$
$\displaystyle<\frac{m_{0}}{\mathcal{S}_{\gamma}(\mathcal{X})},$
$\displaystyle s$
$\displaystyle<\frac{a}{2m_{0}}\Big{[}C\frac{\sigma(\sqrt{D}+\sqrt{d})}{\sqrt{\lambda}}\frac{1}{\sqrt{2s\mathcal{S}_{\gamma}(\mathcal{X})m_{0}-(s\mathcal{S}_{\gamma}(\mathcal{X}))^{2}}}+\frac{c(\mathcal{X})}{\mathcal{S}_{\gamma}(\mathcal{X})}\Big{]}^{-1}.$

The second condition can be satisfied for $T$ greater than a constant
$C^{\prime}$ with respect to $a$. Indeed,

$\Big{[}\frac{1}{\sqrt{sm_{0}}}C\frac{\sigma(\sqrt{D}+\sqrt{d})}{\sqrt{\lambda}}\frac{1}{\sqrt{2\mathcal{S}_{\gamma}(\mathcal{X})-s(\mathcal{S}_{\gamma}(\mathcal{X}))^{2}/m_{0}}}+\frac{c(\mathcal{X})}{\mathcal{S}_{\gamma}(\mathcal{X})}\Big{]}\\\
\leq\Big{[}\frac{1}{\sqrt{sm_{0}}}C\frac{\sigma(\sqrt{D}+\sqrt{d})}{\sqrt{\lambda}}\sqrt{\frac{1}{\mathcal{S}_{\gamma}(\mathcal{X})}}+\frac{c(\mathcal{X})}{\mathcal{S}_{\gamma}(\mathcal{X})}\Big{]}=O(1/\sqrt{sm_{0}}).$

Thus, to satisfy the second condition, we need

$s=\frac{Ca}{T^{\nu}}<C^{\prime}\frac{a}{2m_{0}}\sqrt{m_{0}}\sqrt{\frac{Ca}{T^{\nu}}},$

or

$T=\Omega\Big{(}\Big{[}\frac{m_{0}}{a}\Big{]}^{1/\nu}\Big{)}.$ (B.6)

Choosing $T$ to be sufficiently large so that it also satisfies

$T>\Big{[}\mathcal{C}\frac{m_{0}\log(a/2d_{r}^{2}(\bm{V}_{\star}^{\top}\bm{V}_{0}))}{a\mathcal{S}_{\gamma}(\mathcal{X})}\Big{]}^{1/(1-\nu)}>\frac{\log(a/2d_{r}^{2}(\bm{V}_{\star}^{\top}\bm{V}_{0}))}{\log(1-s\mathcal{S}_{\gamma}(\mathcal{X})/m_{0})},$
(B.7)

then, with probability at least $1-2\lambda$,

$d_{r}^{2}(\bm{V}_{\star},\bm{V}_{T})<a.$ (B.8)

Notice that the constraints on $T$ given by (B.5), (B.6), and (B.7) depend
only on $a$ and $d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0})$ the ratio
$a/d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0})$.

∎

####  B.6.3 Proof of Theorem 3

Since $\mathcal{S}(\mathcal{X}^{k})$ is bounded between
$[-\max_{i}\|\bm{x}_{i}\|^{2},\max_{i}\|\bm{x}_{i}\|^{2}]$, we also have the
uniform bound

$\mathbb{E}[\mathcal{S}(\mathcal{X}^{k})-\mathcal{S}_{\mathbb{E}}]^{2}\leq\max_{i}\|\bm{x}_{i}\|^{2},$
(B.9)

although in general we expect this to be much smaller. In particular, if the
data is spherized, then this is bounded by 1.

###### Proof.

The proof of the theorem follows from the same reasoning as dp-GGD after
splitting the sequence of errors as

$\displaystyle d_{r}^{2}(\bm{V}_{\star},\bm{V}_{k+1})$ $\displaystyle\leq
d_{r}^{2}(\bm{V}_{\star}^{\top}\bm{V}_{k})+s\Big{(}\bm{\beta}_{1}^{\top}\bm{V}_{\star}^{\top}((\mathsf{grad}F(\bm{V}_{k};\mathcal{X}))\bm{\beta}_{2}\Big{)}+c(\mathcal{X})s^{2}$
$\displaystyle\leq\Big{(}1-\frac{s\mathcal{S}_{\mathbb{E}}}{d_{r}(\bm{V}_{\star}^{\top}\bm{V}_{k})}\Big{)}d_{r}^{2}(\bm{V}_{\star},\bm{V}_{k})+(\mathcal{S}_{\gamma}(\mathcal{X}^{k})-\mathcal{S}_{\gamma,\mathbb{E}})sd_{r}(\bm{V}_{\star},\bm{V}_{k})+c(\mathcal{X})s^{2}$
$\displaystyle\leq\Big{(}1-\frac{s\mathcal{S}_{\mathbb{E}}}{m_{T}}\Big{)}d_{r}^{2}(\bm{V}_{\star},\bm{V}_{k})+(\mathcal{S}_{\gamma}(\mathcal{X}^{k})-\mathcal{S}_{\gamma,\mathbb{E}})sd_{r}(\bm{V}_{\star},\bm{V}_{k})+c(\mathcal{X})s^{2},$

and then controlling
$(\mathcal{S}_{\gamma}(\mathcal{X}^{k})-\mathcal{S}_{\gamma,\mathbb{E}})sd_{r}(\bm{V}_{\star},\bm{V}_{k})$.
Here, again, $m_{T}=\max_{j=1,\dots,T}d_{r}(\bm{V}_{\star},\bm{V}_{j})$.

Bounding $m_{T}$: As before, be begin by bounding $m_{T}$ by first looking at
the looser bound

$d_{r}^{2}(\bm{V}_{\star},\bm{V}_{k+1})\leq\Big{(}1-s\mathcal{S}_{\mathbb{E}}\Big{)}d_{r}^{2}(\bm{V}_{\star},\bm{V}_{k})+(\mathcal{S}_{\gamma}(\mathcal{X}^{k})-\mathcal{S}_{\gamma,\mathbb{E}})sd_{r}(\bm{V}_{\star},\bm{V}_{k})+c(\mathcal{X})s^{2}.$

Telescoping yields

$d_{r}^{2}(\bm{V}_{\star},\bm{V}_{k})\leq\Big{(}1-s\mathcal{S}_{\mathbb{E}}\Big{)}^{k}d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0})\\\
+\sum_{j=1}^{k}\Big{(}1-s\mathcal{S}_{\mathbb{E}}\Big{)}^{k-j}(\mathcal{S}_{\gamma}(\mathcal{X}^{j})-\mathcal{S}_{\gamma,\mathbb{E}})sd_{r}(\bm{V}_{\star},\bm{V}_{j})+\sum_{j=1}^{k}\Big{(}1-s\mathcal{S}_{\mathbb{E}}\Big{)}^{k-j}c(\mathcal{X})s^{2}.$

The last term is bounded by

$\sum_{j=1}^{k}\Big{(}1-s\mathcal{S}_{\mathbb{E}}\Big{)}^{k-j}c(\mathcal{X})s^{2}\leq\frac{sc(\mathcal{X})}{\mathcal{S}_{\gamma,\mathbb{E}}}.$

The other term can be bounded by Doob’s maximal inequality

$\Pr\Big{(}\max_{k=1,\dots,T}\Big{|}\sum_{j=1}^{k}\Big{(}1-s\mathcal{S}_{\mathbb{E}}\Big{)}^{k-j}(\mathcal{S}_{\gamma}(\mathcal{X}^{j})-\mathcal{S}_{\gamma,\mathbb{E}})sd_{r}(\bm{V}_{\star},\bm{V}_{j})\Big{|}>\epsilon\Big{)}\\\
\leq\Pr\Big{(}\max_{k=1,\dots,T}\Big{|}\sum_{j=1}^{k}(\mathcal{S}_{\gamma}(\mathcal{X}^{j})-\mathcal{S}_{\gamma,\mathbb{E}})s\Big{|}>\epsilon\Big{)}\leq\frac{s^{2}\sum_{j=1}^{T}\mathbb{E}(\mathcal{S}_{\gamma}(\mathcal{X}^{j})-\mathcal{S}_{\gamma,\mathbb{E}})}{\epsilon^{2}}\leq\frac{s^{2}T}{\epsilon^{2}}.$

Setting $\epsilon=\frac{s\sqrt{T}}{\sqrt{\lambda}}$ yields

$\Pr\Big{(}\max_{k=1,\dots,T}\Big{|}\sum_{j=1}^{k}\Big{(}1-s\mathcal{S}_{\mathbb{E}}\Big{)}^{k-j}(\mathcal{S}_{\gamma}(\mathcal{X}^{j})-\mathcal{S}_{\gamma,\mathbb{E}})sd_{r}(\bm{V}_{\star},\bm{V}_{j})\Big{|}>\frac{s\sqrt{T}}{\sqrt{\lambda}}\Big{)}\leq\lambda.$

Thus, if

$\frac{s\sqrt{T}}{\sqrt{\lambda}}<\frac{d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0})}{2}$

and

$\frac{s^{2}Tc(\mathcal{X})}{\mathcal{S}_{\gamma,\mathbb{E}}}<\frac{d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0})}{2},$

then $m_{T}<2d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0})$. In particular, for
$s=c_{1}aT^{-\nu}$, these are satisfied if

$T>\max\Big{(}\Big{[}\frac{2c_{1}a^{2}c(\mathcal{X})}{\mathcal{S}_{\gamma,\mathbb{E}}d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0})}\Big{]}^{1/(2\nu-1)},\Big{[}\frac{2a}{\sqrt{\lambda}d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0})}\Big{]}^{2/(2\nu-1)}\Big{)}=:\mathcal{F}_{2}(a/d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0}),\lambda).$
(B.10)

Bounding the $T$th iterate: Let $m_{0}=2d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0})$.
From here the proof is straightforward: the first term geometrically
decreases. The second can be bounded with Doob’s maximal inequality with high
probability and uses the fact that

$\sum_{j}(1-s\mathcal{S}_{\gamma,\mathbb{E}}/m_{0})^{2(T-j)}\leq\frac{1}{1-(1-s\mathcal{S}_{\gamma,\mathbb{E}}/m_{0})^{2}}=\frac{m_{0}^{2}}{2s\mathcal{S}_{\gamma,\mathbb{E}}m_{0}-(s\mathcal{S}_{\gamma,\mathbb{E}})^{2}},$

which is independent of $T$. More specifically, Doob’s maximal inequality
yields

$\displaystyle\Pr\Big{(}$ $\displaystyle\inf_{1\leq k\leq
T}\Big{|}\sum_{j=1}^{k}s(1-s\mathcal{S}_{\gamma}(\mathcal{X})/m_{0})^{T-j}(\mathcal{S}_{\gamma}(\mathcal{X}^{j})-\mathcal{S}_{\gamma,\mathbb{E}})d_{r}(\bm{V}_{\star},\bm{V}_{j})\Big{|}>\epsilon\Big{)}$
$\displaystyle\leq\Pr\Big{(}\inf_{1\leq k\leq
T}\Big{|}\sum_{j=1}^{k}s(1-s\mathcal{S}_{\gamma}(\mathcal{X})/m_{0})^{T-j}(\mathcal{S}_{\gamma}(\mathcal{X}^{j})-\mathcal{S}_{\gamma,\mathbb{E}})\Big{|}>\epsilon\Big{)}$
$\displaystyle\leq\frac{\mathbb{E}\Big{(}\sum_{j=1}^{T}s^{2}(1-s\mathcal{S}_{\gamma}(\mathcal{X})/m_{0})^{2(T-j)}\Big{[}(\mathcal{S}_{\gamma}(\mathcal{X}^{j})-\mathcal{S}_{\gamma,\mathbb{E}})\Big{]}^{2}\Big{)}}{\epsilon^{2}}.$

We can upper bound

$\displaystyle\mathbb{E}\Big{[}(\mathcal{S}_{\gamma}(\mathcal{X}^{j})-\mathcal{S}_{\gamma,\mathbb{E}})\Big{]}^{2}$
$\displaystyle\leq 1.$

In any case, this implies that

$\Pr\Big{(}\inf_{1\leq k\leq
T}\Big{|}\sum_{j=1}^{k}s(1-s\mathcal{S}_{\gamma}(\mathcal{X})/m_{0})^{T-j}(\mathcal{S}_{\gamma}(\mathcal{X}^{j})-\mathcal{S}_{\gamma,\mathbb{E}})d_{r}(\bm{V}_{\star},\bm{V}_{j})\Big{|}>\epsilon\Big{)}\\\
\leq\frac{m_{0}^{2}}{2s\mathcal{S}_{\gamma,\mathbb{E}}m_{0}-(s\mathcal{S}_{\gamma,\mathbb{E}})^{2}}\frac{s^{2}}{\epsilon^{2}},$

or

$\Pr\Big{(}\inf_{1\leq k\leq
T}\Big{|}\sum_{j=1}^{k}s(1-s\mathcal{S}_{\gamma}(\mathcal{X})/m_{0})^{T-j}(\mathcal{S}_{\gamma}(\mathcal{X}^{j})-\mathcal{S}_{\gamma,\mathbb{E}})d_{r}(\bm{V}_{\star},\bm{V}_{j})\Big{|}\\\
>s\frac{m_{0}}{\sqrt{2s\mathcal{S}_{\gamma,\mathbb{E}}m_{0}-(s\mathcal{S}_{\gamma,\mathbb{E}})^{2}}\sqrt{\lambda}}\Big{)}\leq\lambda.$

Putting these together, we find with probability at least $1-2\lambda$,

$d_{r}^{2}(\bm{V}_{\star},\bm{V}_{T})\leq(1-s\mathcal{S}_{\gamma}(\mathcal{X})/m_{0})^{T}d_{r}^{2}(\bm{V}_{\star}^{\top}\bm{V}_{0})+s\Big{[}\frac{m_{0}}{\sqrt{2s\mathcal{S}_{\gamma,\mathbb{E}}m_{0}-(s\mathcal{S}_{\gamma,\mathbb{E}})^{2}}\sqrt{\lambda}}+\frac{m_{0}c(\mathcal{X})}{\mathcal{S}_{\gamma,\mathbb{E}}}\Big{]}.$

Now, if $T$ is sufficiently large so that $s=\frac{Ca}{T^{\nu}}$ satisfies

$\displaystyle s$
$\displaystyle<\frac{m_{0}}{\mathcal{S}_{\gamma}(\mathcal{X})}$ $\displaystyle
s$
$\displaystyle<\frac{a}{2m_{0}}\Big{[}\frac{1}{\sqrt{2s\mathcal{S}_{\gamma,\mathbb{E}}m_{0}-(s\mathcal{S}_{\gamma,\mathbb{E}})^{2}}\sqrt{\lambda}}+\frac{c(\mathcal{X})}{\mathcal{S}_{\gamma,\mathbb{E}}}\Big{]}^{-1}.$

To satisfy the second condition, we again need

$T=\Omega\Big{(}\Big{[}\frac{m_{0}}{a}\Big{]}^{1/\nu}\Big{)}.$ (B.11)

Choosing $T$ to be sufficiently large so that it also satisfies

$T>\Big{[}\mathcal{C}\frac{m_{0}\log(a/2d_{r}^{2}(\bm{V}_{\star}^{\top}\bm{V}_{0}))}{a\mathcal{S}_{\gamma,\mathbb{E}}}\Big{]}^{1/(1-\nu)}>\frac{\log(a/2d_{r}^{2}(\bm{V}_{\star}^{\top}\bm{V}_{0}))}{\log(1-s\mathcal{S}_{\gamma,\mathbb{E}}/m_{0})},$
(B.12)

then, with probability at least $1-2\lambda$,

$d_{r}^{2}(\bm{V}_{\star},\bm{V}_{T})<a.$ (B.13)

Notice again that the constraints on $T$ given by (B.5), (B.6), and (B.7)
depend only on $a$ and $d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0})$ the ratio
$a/d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0})$. ∎

Finally, we demonstrate how one might hope to have
$\mathcal{S}_{\gamma,\mathbb{E}}>0$ with a simple model. If the inlier and
outlier percentages are $\alpha_{1}$ and $\alpha_{0}$ and the inlier and
outlier sample covariances are $\bm{\Sigma}_{\mathrm{in}}$ and
$\bm{\Sigma}_{\mathrm{out}}$, respectively, then

$\mathbb{E}\mathcal{S}(\mathcal{X}^{k})\geq\mathbb{E}\Big{(}\gamma\lambda_{r}\Big{(}\frac{1}{m}\bm{X}_{\mathrm{in}}^{k}\bm{X}_{\mathrm{in}}^{k\top}\Big{)}-\lambda_{1}\Big{(}\frac{1}{m}\bm{X}_{\mathrm{out}}^{k}\bm{X}_{\mathrm{out}}^{k\top}\Big{)}\Big{)}=\gamma\alpha_{0}\lambda_{r}(\bm{\Sigma}_{\mathrm{in}})-\alpha_{1}\lambda_{1}(\bm{\Sigma}_{\mathrm{out}})=\mathcal{S}_{\gamma,\mathbb{E}}$
(B.14)

Therefore, one could assume that
$\gamma\alpha_{0}\lambda_{r}(\bm{\Sigma}_{\mathrm{in}})-\alpha_{1}\lambda_{1}(\bm{\Sigma}_{\mathrm{out}})>0$.

####  B.6.4 Proof of Theorem 4

###### Proof.

Combining the results of the previous two theorems yields the result by
simultaneously controlling both martingales. Notice that $\mathcal{F}_{3}$
will be defined similarly to $\mathcal{F}_{1}$ and $\mathcal{F}_{2}$.

∎

####  B.6.5 Proof of Theorem 5

###### Proof.

Set $s=\frac{c_{1}a}{T^{\nu}}$. For $T_{1}$ sufficiently large so that the
conditions within the theorem hold. Then, with probability $1-2\lambda$ (or
$1-4\lambda$ for dp-SGGD), in $T_{1}$ iterations,
$d_{r}^{2}(\bm{V}_{\star},\bm{V}_{T_{1}})<a$.

Now suppose that we restart with $s^{\prime}=s/2$ and
$\bm{V}_{0}=\bm{V}_{T_{1}}$. Notice that this is equivalent to taking
$a^{\prime}=a/2$ and starting distance
$d_{r}^{2}(\bm{V}_{\star},\bm{V}_{0}^{\prime})=d_{r}^{2}(\bm{V}_{\star},\bm{V}_{T})<a$.
In particular, this takes $T^{\prime}>\mathcal{F}_{\bullet}(1/2)$ iterations
to reach $\bm{V}_{T+T^{\prime}}$ such that
$d_{r}^{2}(\bm{V}_{\star},\bm{V}_{T+T^{\prime}})<a/2$. Repeating this
procedure for $r$ restarts every $\mathcal{F}_{\bullet}(1/2)$ iterations
yields the desired result.

∎

This section gives additional plots demonstrating the performance of the
various differentially private methods we discuss.

Figure 4: $\delta$ and $\epsilon$ phase transition plot. Each square
represents $\log_{10}$ final error.

First, we give a phase transition plot with respect to $\epsilon$ and
$\delta$. Figure 4 shows that dp-GGD’s transition from small $\epsilon$ and
$\delta$ is more abrupt than that of dp-GD-REAP and dp-MD-REAP. The data
parameters are as follows, the inlier dimension $r=2$, total dimension $D=20$,
number of points $N=2,000$ and inlier ratio $=0.5$. The algorithms’ parameters
are as follows, total number of iterations of each algorithm are the same as
the number of data points $N$, the step size for the four dp-REAP algorithms
to be $\eta_{k}=8/\sqrt{k}$, the step size for dp-GGD and dp-SGGD is
$\eta_{k}=1/2^{\lfloor k/50\rfloor}$. The experiment is repeated 50 times and
the medium error is plotted.

Next, we plot the ratio of successful attempts to converge to tolerance
$10^{-2}$ in 50 repetitions of each algorithm as a function of inlier ratio
and batch size in Figure 5. (plot inlier ratio vs N) The data parameters are
as follows, number of points $N=2,000$, total dimension $D=20$, and inlier
dimension $r=2$. The algorithms’ parameters are as follows, total number of
iterations of each algorithm are the same as the number of points $N$, the
step size for the four dp-REAP algorithms to be $\eta_{k}=8/\sqrt{k}$, the
step size for dp-GGD and dp-SGGD is $\eta_{k}=1/2^{\lfloor k/50\rfloor}$. The
plot shows that dp-SGD-REAP converges in the regime where inlier ratio is
$\geq 0.6$, dp-SMD-REAP converges in the regime where inlier ratio is $\geq
0.2$, and dp-SGGD converges almost for all inlier ratios when the batch size
is greater than 2.

Figure 5: Inlier ratio and batch size transition plot. Each square represents
the percentage of time that the algorithm achieves the given tolerance for a
given combination of inlier ratio and batch size.

We plot the percentage of repetitions that algorithms converge to tolerance
$10^{-2}$ in Figure 5.

We also give a phase transition plot of $D$ versus $r$, where the value is the
final $\log_{10}$ error. Figure 6 shows that dp-GD-REAP, dp-MD-REAP and dp-GGD
work well in the regime where both $D$ and $r$ are small. The data parameters
are as follows, number of points $N=2,000$ and inlier ratio $=0.5$. The
algorithms’ parameters are as follows, total number of iterations of each
algorithm are the same as the number of points $N$, the step size for the four
dp-REAP algorithms to be $\eta_{k}=8/\sqrt{k}$, the step size for dp-GGD and
dp-SGGD is $\eta_{k}=1/2^{\lfloor k/50\rfloor}$. The experiment is repeated 50
times and the medium error of all repetitions is plotted.

Figure 6: $D$ and $r$ phase transition plot. Each square represents
$\log_{10}$ final error. The $x$-axis correspond to different values of $D$
and the $y$-axis correspond to different values of $r$.

Next, we give a phase transition plot of $N$ versus $D$ with $\log_{10}$ error
in Figure 7, and time to error = $10^{-2}$ in Figure 8. We plot $\log_{10}$
final error, and time to converge to tolerance (Figure 8), and ratio of failed
attempt to reach tolerance in 50 repetitions (Figure 9) as a function of $N$
and $D$. In the event that none of the repetitions successfully reaches
tolerance, the square shows yellow in Figure 8, this corresponds with number
of failed attempts in Figure 9. Figure 7 shows that the dp-GD-REAP, dp-MD-REAP
and dp-GD work well in the regime where $D$ is small and $N$ is large. The
data parameters are as follows, the inlier dimension $r=2$, and the inlier
ratio $=0.5$. The algorithms’ parameters are as follows, total number of
iterations of each algorithm are the same twice the number of points $2N$, the
step size for the four dp-REAP algorithms to be $\eta_{k}=8/\sqrt{k}$, the
step size for dp-GGD and dp-SGGD is $\eta_{k}=1/2^{\lfloor k/50\rfloor}$. Each
algorithm is run repetitions of 50 times.

Figure 7: $D$ and $N$ phase transition plot. Each square represents
$\log_{10}$ final error.

Figure 8: $D$ and $N$ phase transition plot. Each square represents the time
to reach tolerance $10^{-2}$. If in all repetitions the algorithms failed to
converge, the square is shown in yellow (I imputed a large number, 2 in this
case).

Figure 9: $D$ and $N$ phase transition plot. Each square represents the
percentage of repetitions each algorithm fails to reach tolerance $10^{-2}$.