arXiv:1507.00056v2 [cs.DS] 25 Nov 2015

Private Approximations of the 2nd-Moment Matrix Using Existing
Techniques in Linear Regression
Or Sheffet
Center for Research on Computation and Society
Harvard University
Cambridge, MA
osheffet@seas.harvard.edu
November 26, 2015
Abstract
We introduce three differentially-private algorithms that approximate the 2nd-moment matrix of
the data. These algorithm, which in contrast to existing algorithms output positive-definite matrices,
correspond to existing techniques in linear regression literature. Specifically, we discuss the following
three techniques. (i) For Ridge Regression, we propose setting the regularization coefficient so that by
approximating the solution using Johnson-Lindenstrauss transform we preserve privacy. (ii) We show
that adding a small batch of random samples to our data preserves differential privacy. (iii) We show that
sampling the 2nd-moment matrix from a Bayesian posterior inverse-Wishart distribution is differentially
private provided the prior is set correctly. We also evaluate our techniques experimentally and compare
them to the existing “Analyze Gauss” algorithm of Dwork et al [DTTZ14].

1

Introduction

Differentially private algorithms [DMNS06, DKM+ 06] are data analysis algorithms that give a strong guarantee of privacy, roughly stated as: by adding to or removing from the data a single datapoint we do not
significantly change the probability of any outcome of the algorithm. The focus of this paper is on differentially private approximations of the 2nd-moment matrix of the data — given a dataset D ∈ Rn×d , its
2nd-moment matrix (also referred to as the Gram matrix of data or the scatter matrix if the mean of D is 0)
is the matrix DT D — and the uses of such approximations in linear regression. Indeed, since the 2nd-moment
matrix of the data plays a major role in many data-analysis techniques, we already have differentially private
algorithms that approximate the 2nd-moment matrix [DTTZ14] for the purpose of approximating the PCA,
techniques for approximating the rank-k PCA of the data directly [HR12, Har13, KT13], or differentially
private algorithms for linear regressions [CMS11, KST12, TS13, BST14].
However, existing techniques for differentially private linear regression suffer from the drawback that they
approximate a single regression. That is, they assume that each datapoint is composed of a vector of features
x and a label y and find the best linear combination of the features that predicts y. Yet, given a dataset
D with d attributes we are free to pick any single attribute as a label, and any subset of the remaining
attributes as features. Therefore, a database with d attributes yields exp(d) potential linear regression
problems; and running these algorithms for each linear regression problem separately simply introduces far
too much random noise.1
In contrast, the differentially private techniques that approximate the 2nd-moment matrix of the data,
such as the Analyze Gauss paper of Dwork et al [DTTZ14], allow us to run as many regressions on the data
1 Indeed, Ullman [Ull15] have devised a solution to this problem, but this solution works in the more-cumbersome online
model and requires exponential running-time for the curator; whereas our techniques follow the more efficient offline approach.

1

as we want. Yet, to the best of our knowledge, they have never been analyzed for the purpose of linear
regression. Furthermore, the Analyze Gauss algorithm suffers from the drawback that it does not necessarily
output a positive-definite matrix. This, as discussed in [XKI11] and as we show in our experiments, can be
very detrimental — even if we do project the output back onto the set of positive definite matrices. And
though the focus of this work is on linear regression, one can postulate additional reasons why releasing a
positive definite matrix is of importance, such as using the output as a kernel matrix or doing statistical
inference on top of the linear regression.
Our Contribution. In this work, we give three differentially private techniques for approximating the
2nd-moment matrix of the data that output a positive-definite matrix. We analyze their utility, both
theoretically and empirically, and more importantly — show how they correspond to existing techniques in
linear regression. And so we contribute to an increasing line of works [BBDS12, VZ15, WFS15] that shows
that differential privacy may rise from existing techniques, provided parameters are set properly. We also
compare our algorithms to the existing Analyze Gauss technique.
(Some notation before we introduce our techniques. We assume the data is a matrix A ∈ Rn×d with n
sample points in d dimensions. For the ease of exposition, we focus on a single regression problem, given by
A = [X; y] — i.e., the label is the d-th column and the features are the remaining p = d − 1 columns. We
use σmin (A) to denote the least singular value of A.)
1. The Johnson-Lindenstrauss Transform and Ridge Regression. Blocki et al [BBDS12] have shown
that projecting the data using a Gaussian Johnson-Lindenstrauss transform preserves privacy if σmin (A)
is sufficiently large and it has been applied for linear regression [Upa14]. Our first result improves on the
analysis of Blocki et al and uses a smaller bound on σmin (A) (shaving off a factor of log(r) with r denoting
the number of rows in the JL transform). This result implies that when σmin (A) is large we can project
the data using the JL-transform and output the 2nd-moment matrix of the projected data and preserve
privacy. Furthermore, it is also known [Sar06] that the JL-transform gives a good approximation for linear
regression problems. However, this is somewhat contradictory to our intuition: for datasets where y is
well approximated by a linear combination of X, the least singular value should be small (as A’s stretch
along the direction (β, −1) T is small). That is why we artificially increase the singular values of A by
appending it with a matrix w · Id×d . It turns out that this corresponds to approximating the solution of the
Ridge regression problem [Tik63,
HK70], the linear regression problem with l2 -regularization — the problem
P
of finding β R = arg minβ i kyi − β · xi k2 + w2 kβk2 . Literature suggests many approaches [HTF09] to
determining the penalty coefficient w2 , approaches that are based on the data itself and on minimizing risk.
Here we give a fundamentally different approach — set w as to preserve (, δ)-differential privacy. Details,
utility analysis and experiments regarding this approach appear in Section 3.
2. Additive Wishart noise. Whereas the Analyze Gauss algorithm adds Gaussian noise to AT A, here we
show that we can sample a positive definite matrix W from a suitably chosen Wishart distribution Wd (V, k),
and output AT A + W . This in turn corresponds to appending A with k i.i.d samples from a multivariate
Gaussian N (0d , V ). One is able to view this too as an extension of Ridge regression, where instead of
appending A with d fixed examples, we append A with k ≈ d + O(1/2 ) random examples.2 Note, as
opposed to Analyze Gauss [DTTZ14], where the noise has 0-mean, here the expected value of the noise
is kV . This yields a useful way of post-processing the output: AT A + W − kV . Details, theorems and
experiments with additive Wishart noise appear in Section 4.
3. Sampling from an inverse-Wishart distribution. The Bayesian approach for estimating the 2ndmoment matrix of the data assumes that the n sample points are sampled i.i.d from some N (0d , V ) for
some unknown V , where we have a prior distribution on V . Each sample point causes us to update our
belief on V which results in a posterior distribution on V . Though often one just outputs the MAP of
the posterior belief (the mean of the posterior distribution), it is also common to output a sample drawn
randomly from the posterior distribution. We show that if one uses the inverse-Wishart distribution as
a prior (which is common, as the inverse-Wishart distribution is a conjugate prior), then sampling from
2 Though it is also tempting to think of this technique as running Bayesian regression with random prior, this analogy does
not fully carry through as we discuss later.

2

the posterior is (, δ)-diffrentially private, provided the prior is spread enough. This gives rise to our third
approach of approximating AT A — sampling from a suitable inverse Wishart distribution. We comment
that the idea that existing techniques in Bayesian analysis, and specifically sampling from the posterior
distribution, are differentially-private on their own was originally introduced in the beautiful and elegant
work of Vadhan and Zheng [VZ15]. But whereas their work focuses on estimating the mean of the sample,
we focus on estimating the variances/2nd-moment. Details, theorems and experiments on sampling from the
inverse-Wishart distribution appear in Section 5.
Finally, in Section 6 we compare our algorithms to the Analyze Gauss algorithm. We show that in the
simple case where the data is devised by p independent features concatenated with a single linear combination
of the features, the Analyze Gauss algorithm, which introduces the least noise out of all algorithms, is clearly
the best algorithm once n is sufficiently large. However, when the data contains multiple such regressions
and therefore has small singular values, the situation is far from being clear cut, and indeed, unless n is
extremely large, our algorithms achieve smaller errors than the Analyze Gauss baseline. We comment that
our experiments should be viewed solely as a proof-of-conecpt. They are only preliminary, and much more
experimentation is needed to fully evaluate the benefits of the various algorithms.
Our proof technique. Before continuing to preliminaries and the formal details of our algorithms, we
give an overview of the proof technique. (All of the proofs are deferred to Appendix B.) To prove that each
algorithm preserves (, δ)-differential privacy we state and prove 3 corresponding theorems, whose proofs
follow the same high-level approach. As mentioned above, one theorem improves on a theorem of Blocki et
al [BBDS12], who were the first to show that the JL-transform is differentially private. Blocki et al observed
that by projecting the data using a (r × n)-matrix of i.i.d normal Gaussians, we effectively repeat the same
one-dimensional projection r independent times. So they proved that each one-dimensional projection is
(, δ)-differentially private, and to show the entire projection preserves privacy
they used the off-the-shelf
√
r
log(r)).
In order to derive a
composition of Dwork et al [DRV10],
getting
a
bound
that
depends
on
O(
√
bound depending only on O( r), we do not use the composition theorem of [DRV10] but rather study the
specific r-fold composition of the projection. As a result, we cannot follow the approach of Blocki et al.
To show that a one-dimensional projection is (, δ)-differentially private, Blocki et al compared the PDFs
of two multivariate Gaussians. The PDF of a multivariate Gaussian is given by the multiplication of two
terms: the first depends on the determinant of the variance, and the second depends on some exponent (see
exact definition in Section 2). Blocki et al compared the ratio of each of the terms and showed that w.h.p
each term’s ratio is bounded by e/2 . Unfortunately, following the same approach of Blocki et al yields a
bound of er/2 for each of the terms and an overall bound that depends on O(r). Instead, we observe that
the contributions of the determinant term and the exponent term to the ratio of the PDFs are of opposite
signs. So we use the Matrix Determinant Lemma and the Sherman-Morrison Lemma (see Theorem A.4)
to combine both terms into a single exponent term, and bound its size using the Johnson-Lindenstrauss
transform (or rather, tight bounds on the χ2 -distribution). The main lemma we use in our analysis is
detailed in Lemma A.1. This lemma, in addition to giving tight bounds for the Gaussian JL-transform
(mimicking the approach of Dasgupta and Gupta [DG03]), also gives a result that might be of independent
interest. The standard JL lemma shows that for a (r × d)-matrix R of i.i.d normal Gaussians and any fixed
vector v it holds w.h.p that v T v ∈ (1 ± η)v T ( 1r RT R)v provided r = O(η −2 ). In Lemma A.1 we also show
1
that for any fixed v we have w.h.p. that v T v ∈ (1 ± η)v T ( r−d
RT R)−1 v provided r = d + O(η −2 ). 3

2

Preliminaries and Notation

Notation. Throughout this paper, we use lower-case letters to denote scalars; bold characters to denote
vectors; and UPPER-case letters to denote matrices. The l-dimensional all zero vector is denoted 0l , and
the (l × m)-matrix of all zeros is denoted 0l×m . The l-dimensional identity matrix is denoted Il×l . For two
3 To the best of our knowledge, for a general JLT, this is known to hold only when r = O(d · η −2 ) and the transform preserves
the lengths of all vectors in the Rd space, see [Sar06] Corollary 11.

3

matrices M, N with the same number of row we use [M ; N ] to denote the concatenation of M and N . We
use , δ to denote the privacy parameters. P
For a given matrix, kM k denotes the spectral norm (= σmax (M ))
2 1/2
and kM kF denotes the Frobenious norm ( j,k Mj,k
) ; and use σmax (M ) and σmin (M ) to denote its largest
and smallest singular value resp.
The Gaussian Distribution and Related Distributions. We denote by Lap(σ)
the Laplace dis
tribution whose mean is 0 and variance is 2σ 2 . A univariate Gaussian N µ, σ 2 denotes the Gaussian
distribution whose
mean is µ and variance σ 2 . Standard concentration bounds on Gaussians give that
p
Pr[x > µ + σ ln(1/ν)] < ν. A multivariate Gaussian N (µ, Σ) for some positive semi-definite Σ denotes
the multivariate Gaussian distribution where the mean of the j-th coordinate is the µj and the co-variance
between coordinates j and k is Σj,k . The PDF of such Gaussian is defined only on the subspace colspan(Σ),


−1/2

T
˜
where for every x ∈ colspan(Σ) we have PDF(x) = (2π)rank(Σ) · det(Σ)
exp − 12 (x − µ) Σ† (x − µ)
˜
and det(Σ)
is the multiplication of all non-zero singular values of Σ. We will repeatedly
re use the rules

2
2 2
garding linear operations on Gaussians. That is, for
any
scalar
c,
it
holds
that
cN
µ,
σ
=
N
c
·
µ,
c
σ
.


For any matrix C it holds that C · N (µ, Σ) = N Cµ, CΣC

T

.

χ2k -distribution,

The
where k is referred to as the degrees of freedom of the distribution, is the distribution
over the l2 -norm of the sum of k independent normal Gaussians. That is, given X1 , . . . , Xk ∼ N (0, 1)
def

it holds that ζ = (X1 , X2 , . . . , Xk ) ∼ N (0k , Ik×k ), and kζk2 ∼ χ2k . Standard tail bounds on the χ2 2
√
p
distribution give that for any ν ∈ (0, 21 ) we have Prx∼χ2k [x ∈
k ± 2 ln(2/ν) ] ≥ 1 − ν. (We present
them in Section A for completeness.) The Wishart-distribution Wd (V, m) is the multivariate extension of
the χ2 -distribution. It describes the scatter matrix of a sample of m i.i.d samples from a multivariate
Gaussian N (0d , V ) and so the support of the distribution is on positive definite matrices. For m > d − 1 we
m−d−1
m
have that PDFWd (V,m) (X) ∝ det(V )− 2 det(X) 2 exp(− 12 tr(V −1 X). The inverse-Wishart distribution
Wd−1 (V, m) describes the distribution over positive definite matrices whose inverse is sampled from the
Wishart distribution using the inverse of V ; i.e. X ∼ Wd−1 (V, m) iff X −1 ∼ Wd (V −1 , m). For m > d − 1 it
m+d+1
m
holds that PDFW −1 (V,m) (X) ∝ det(V ) 2 det(X)− 2 exp(− 12 tr(V X −1 )).
d

Differential Privacy. In this work, we deal with input of the form of a (n × d)-matrix with each row
bounded by a l2 -norm of B. Converting A into a linear regression problem, we denote A as the concatenation
of the (n × p)-matrix X with the vector y ∈ Rn (A = [X; y]) where p = d − 1. This implies we are tying to
predict y as a linear combination of the columns of X. Two matrices A and A0 are called neighbors if they
differ on a single row.
Definition 2.1 ([DMNS06, DKM+ 06]). An algorithm ALG which maps (n × d)-matrices into some range
R is (, δ)-differential privacy if for all pairs of neighboring inputs A and A0 and all subsets S ⊂ R it holds
that Pr[ALG(A) ∈ S] ≤ e Pr[ALG(A0 ) ∈ S] + δ. When δ = 0 we say the algorithm is -differentially private.
It was shown in [DMNS06] that for any f where kf (A) − f (A0 )k1 ≤ ∆ then the algorithm that adds
Laplace noise Lap( ∆ ) to f (A) is -differential privacy.
It wasshown in [DKM+ 06] that for any f where

2
kf (A)−f (A0 )k2 ≤ ∆ then adding Laplace noise N 0, 2∆ ln(2/δ)
to f (A) is (, δ)-differential privacy. This is

precisely the algorithm of Dwork et al in their “Analyze Gauss” paper [DTTZ14]. They observed that in our
setting, for the function f (A) = AT A we have that kf (A) − f (A0 )k2F = B 4 . And so they add i.i.d Gaussian
noise to each coordinate of AT A (forcing the noise to be symmetric, as AT A is symmetric). We therefore
refer to this benchmark as the Analyze Gauss algorithm. In addition, it is known that the composition of
two algorithms, each of which is (, δ)-differentially private, yields an algorithm which is (2, 2δ)-differentially
private.

4

3

Ridge Regression — Set the Regularization Coefficient to Preserve Privacy

The standard problem of linear regression, finding βb = arg minβ kXβ − yk2 , relies on the fact that X is of
full-rank. This clearly isn’t always the case, and X T X may be singular or close to singular. To that end,
as well as for the purpose of preventing over-fitting, regularization is introduced. One way to regularize the
linear regression problem is to introduce a l2 -penalty term: finding β R = arg minβ kXβ − yk2 + w2 kβk2 .
This is known as the Ridge regression problem, introduce by [Tik63, HK70] in the 60s and 70s. Ridge
regression has a closed form solution: β R = (X T X + w2 Ip×p )X T y. The problem of setting w has been
well-studied [HTF09] where existing techniques are data-driven, often proposing to set w as to minimize the
risk of β R . Here, we propose a fundamentally different approach to the problem of setting w: set it so that
we can satisfy (, δ)-differential privacy (via the Johnson-Lindenstrauss transform).
Observe, the Ridge regression problem can be written as: minimizekXβ − yk2 + kwIp×p β − 0p k2 . So,
denote X 0 are the ((n + p) × p)-matrix which we get by concatenating X and wIp×p , and denote y 0 as
the concatenation of y with p zeros. Then β R = arg min kX 0 β − y 0 k2 . Since p = d − 1 and we denote
A = [X; y], we can in fact set A0 as the concatenation of A with the d-dimensional matrix wId×d , and we

 2
β
def
= kX 0 β − y 0 k2 + w2 . Hence β R = arg min f (β). Hence, an approximation
have that f (β) = A0
−1
of A0T A0 yields an approximation of the Ridge regression problem. One way to approximate A0T A0 is via the
Johnson-Lindenstrauss transform, which is known to be differentially private if all the singular values of the
given input are sufficiently large [BBDS12]. And that is precisely why we use A0 — all the singular values
of A0T A0 are greater by w2 than the singular values ofAT A, and in particular are always ≥ w2 . Therefore,
applying the JLT to A0 gives an approximation of A0T A0 , and furthermore, due to the work of Sarlos [Sar06]
the JLT also approximates the linear regression. The following theorem improves on the original theorem of
Blocki et al [BBDS12].
Theorem 3.1. Fix  > 0 and δ ∈ (0, 1e ). Fix B > 0. Fix a positive integer r and let w be such that
q

w2 = 4B 2
2r ln( 4δ ) + ln( 4δ ) /. Let A be a (n × d)-matrix with d < r and where each row of A has
bounded L2 -norm of B. Given that σmin (A) ≥ w, the algorithm that picks a (r × n)-matrix R whose entries
are i.i.d samples from a normal distribution N (0, 1) and publishes R · A is (, δ)-differentially private.
This gives rise to our first algorithm. Algorithm 1 gets as input the parameter r — the number of rows
in our JLT, and chooses the appropriate regularization coefficient w. Based on Theorem 3.1 and abovementioned discussion, it is clear that Algorithm 1 is (, δ)-differentially private. Furthermore, based on the
work of Sarlos, we can also argue the following.
Theorem 3.2. [[Sar06], Theorem 12] Fix any η > 0 and ν ∈ (0, 21 ). Apply Algorithm 1 with r =
O(d log(d) ln(1/ν)/η 2 ). Then, w.p ≥ 1 − ν it holds that kβ R − βeR k ≤ √ 2 η
f (β R ).
T
w +σmin (A A)

b 2 ] (see [DFKU13]) can be used together with
Existing results about the expected distance E[kβ R − βk
R
2
b .
Theorem 3.2 to give a bound on kβe − βk
In addition to Algorithm 1, we can use part of the privacy budget to look at the least singular-value of
AT A. If it happens to be the case that σmin (AT A) is large, then we can adjust w by decreasing it by the
appropriate factor. In fact, one can completely invert the algorithm and, in case σmin (AT A) is really large,
not only set the regularization coefficient to be any arbitrary non-negative number, but also determine r
based on Thm 3.1. Details appear in Algorithm 2.
To measure the effect of regularization we ran the following experiment. (Since the same experimental
setting is used in the following sections we describe it here lengthly, and refer to it in later sections.)

5

Input: A matrix A ∈ Rn×d and a bound B > 0 on the l2 -norm of any row in A.
Privacy parameters: , δ > 0.
Parameter
r indicating the number of rows in the resulting matrix.
r

q
2r ln( 4δ ) + ln( 4δ ) /.
Set w = 4B 2
Set A0 as the concatenation of A with wId×d .
Sample a r × (n + d)-matrix R whose entries are i.i.d samples from a normal Gaussian.
return M = 1r (RA0 )T (RA0 ) and the approximation βeR = arg minβd =−1 β T M β.
Algorithm 1: Approximating Ridge Regression while Preserving Privacy
Input: A matrix A ∈ Rn×d and a bound B > 0 on the l2 -norm of any row in A.
Privacy parameters: , δ > 0.
Parameter
r0 indicating the minimal number of rows in the resulting matrix.
r
q

8
8B 2
8
2r
ln(
Set w =
)
+
ln(
)
.
0

δ
δ
o
n
2
2
+ Z where Z ∼ Lap( 2B ).
Set s ← max 0, σmin (AT A) − 2B ln(2/δ)

p
Adjust w ← max{0, w2 − s}.
if w > 0 then
Set A0 as the concatenation of A with wId×d .
Sample a r × (n + d)-matrix R whose entries are i.i.d samples from a normal Gaussian.
return M = r10 (RA0 )T (RA0 ), w and the approximation βeR = arg minβd =−1 β T M β.
else

o
n
q
Set r∗ ← max r ∈ Z :

8B 2


2r ln( 8δ ) + ln( 8δ ) ≤ s .

Sample a (r∗ × n)-matrix R whose entries are i.i.d samples from a normal Gaussian.
return M = r1∗ (RA)T (RA), r∗ and the approximation βe = arg minβd =−1 β T M β.
end
Algorithm 2: Approximating Regression (Ridge or standard) while Preserving Privacy.

3.1

The Basic Single-Regression Experiment — Setting

To compare between the various algorithms we introduce and to analyze their utility we ran experiments
testing their performance over data generated from a multivariate Gaussian. The experiments all share the
same common setting, but each experiment studied a different set of estimators. In this section we detail
the common setting, and in the next one we details the specific estimators and results of each experiment
separately.
We pick p = 20 i.i.d. features sampled from a normal Gaussian, and pick some β ∈R [−1, 1]p+1 (the
last coordinate denotes the regression’s intercept), and set y as the linear combination of the features and
the intercept (the all-1 column) plus random noise sampled from N (0, 0.5). Hence our data had dimension
d = p+2 = 22 and the 21-dimensional vector β has l2 of about 3. We vary n to take any of the values in {214 =
4, 096, 215 , 216 , . . . , 225 = 33, 554, 432}. We vary
√  to take any of the values {0.05, 0.1, 0.15, 0.2, 0.25, 0.5}, and
fix4 δ = e−9 , and use the l2 -bound of B = 2.5d. (As preprocessing, each datapoint whose length is > B
is shrunk to have length B.) For each estimator we experimented with, we run it t = 15 times, and report
the mean and standard variation of the 15 experiments. In all experiments we measure the l2 -distance
between the outputted estimator of each algorithm to the true β we used to generate the data. After all,
the algorithms we give are aimed at learning the β that generated the given samples, and so they should
4 We are aware that it is a good standard practice to set δ < 1 since otherwise, sampling from the data is (, δ)-differentially
n
private. However, as we vary n drastically, we aim to keep all other parameters equal.

6

return an estimator close to the true β. We coded all experiments in R and ran the experiments on standard
laptop.

3.2

Experiment on Ridge Regression — Measuring the Effect of Regularization

To measure the effect of regularization we ran the following experiment in the setting detailed in Section 3.1.
For each choice of  and n we ran three predictors. The first one is based on Algorithm 2 with r0 = 2d. The
second one is Algorithm 1 where we fed it the parameter r that the first one used, just so all predictors will
be comparable. The last one is the non-private version that projected the data itself, without appending it
with the w · Id×d matrix (again, using the same parameter r as the other two predictors). The results are
given in Figure 1.

Figure 1: (best seen in color) A comparison of the average l2 -error of the JL-based estimators. Algorithm 2
in blue, Algorithm 1 in red, and the non-private version in black. The x-axis is the size of the data in
log-scale.

The results are strikingly similar across all values of . Initially the error of the predictors is very high
(for the value β we used to generate the data, kβk ≈ 2.786, so such levels on noise mean in fact zero utility).
Furthermore, it takes a while until Algorithm 2 (in blue) outperforms the more naı̈ve Algorithm 1 (in red).
(In most experiments, it happens only once n ≥ 219 or n ≥ 220 .) This implies that the privacy-budget
“wasted” on the private estimation of the least singular value of the data actually ends up reducing our
utility but not by a large factor. Towards the largest value of n, Algorithm 2 actually does noticeably better
than Algorithm 1 by a multiplicative factor of ≈ 3.5 to ≈ 11 (for n = 225 when  = 0.1 we have mean accuracy
of 0.0192 vs. 0.0671; when  = 0.5 we have mean accuracy of 0.0058 vs. 0.0639). In all experiments, the
non-private estimator (in black) was clearly the best for all values of n.

7

4

Additive Wishart Noise — Regression with Additional Random
Examples

As discussed in the previous section, Ridge regression can be viewed as regression where in addition to
the sample points given by [X; y] we see d additional datapoints given by wId×d . Our second techniques
follows this approach, only, instead of introducing these d fixed datapoints, we introduce a few more than
d datapoints which are random and independent of the data.5 Formally, we give the details in Algorithm 3
and immediately following — the theorem proving it is (, δ)-differentially private.
Input: A matrix A ∈ Rn×d and a bound B > 0 on the l2 -norm of any row in A.
Privacy parameters: , δ > 0.
Set k ← bd + 14
2 · 2 ln(4/δ)c.

Sample v1 , v2 , . . . , vk i.i.d examples from N 0d , B 2 Id×d .
Pk
return M = AT A + i=1 vi vi T and the approximation βe = arg minβd =−1 β T M β.
Algorithm 3: Additive Wishart Noise Algorithm
Theorem 4.1. Fix  ∈ (0, 1) and δ ∈ (0, 1e ). Fix B > 0. Let A be a (n × d)-matrix where each row of A
has bounded l2 -norm of B. Let N be a matrix sampled from the d-dimensional Wishart distribution with
k-degrees of freedom using the scale matrix B 2 · Id×d (i.e., N ∼ Wd (B 2 · Id×d , k)) for k ≥ bd + 14
2 · 2 ln(4/δ)c.
Then outputting X = AT A + N is (, δ)-differentially private.
Note: Ridge Regression also has a Bayesian interpretation, as introducing a prior on β in regression
problem. It is therefore tempting to argue that Theorem 4.1 implies that solving the regression problem
with a random prior preserves privacy. (I.e., output the MAP of β after setting its prior to a random sample
from the Wishart distribution.) However, this analogy isn’t fully accurate, since our algorithm also adds
random noise to X T y. Indeed, regardless of what prior we use for β, if y = 0n then we always output 0p as
the estimator of β, so one can differentiate between the case that y = 0n and y 6= 0n . We leave the (very
interesting) question of whether Wishart additive random noise can be interpreted as a Bayesian prior for
future work.
We give a bound on the utility of the estimator we get with this technique in Theorem C.3. However,
we are more interested in the utility of this approach after we remove some of the noise we add in this
technique. Note, E[N ] = kB 2 · Id×d , and so it stands to reason that we output AT A + N − kB 2 · Id×d .
Now, when σmin (AT A) is small, we run the risk that some of the eigenvalues of AT A + N are smaller
then kB 2 , causing some of the eigenvalue of AT A + N − kB 2 · Id×d to be negative (which means we no
longer output a PSD). In such a case, Lemma A.3 assures us that w.h.p we can decrease AT A + N by
√
2
p
√
B2
k − ( d + 2 ln(4/δ)) · Id×d and maintain the property that the output is positive definite matrix.
This is the algorithm we set to evaluate empirically.

4.1

Experiment on Additive Wishart Noise

To evaluate the utility of the additive random Wishart noise algorithm we implemented and ran the algorithm
in the same setting as detailed in Section 3.1. For each choice of  and n we ran three predictors. The first one
b The second
is the naı̈ve and non-private linear regression, that uses the data with no additive noise (i.e., β).
one is given by Algorithm 3. The last one is the estimator we get using the output of Algorithm 3 minus
√
2
p
√
either kB 2 · Id×d or B 2
k − ( d + 2 ln(4/δ)) · Id×d (whichever of the two we can use and maintain
positive definiteness). We repeat each experiment t = 15 times, measuring the l2 -distance between the
outputted estimator of each algorithm to the true β we used to generate the data. (This yields randomness
5 Independent

of the data itself, but dependent of its properties. Our noise does depend on the l2 -bound B.

8

in kβb − βk, since every time we re-sample the data.) We report the mean and standard variation of the 15
experiments. The results are given in Figure 2. The results are again consistent across the board — reducing
the noise also reduces the error, and indeed the second estimator is consistently doing better than the naı̈ve
estimator.

Figure 2: (best seen in color) A comparison of the average l2 -error of the Wishart additive noise estimators.
Algorithm 3 in blue, deducting the expected shift from the output of Algorithm 3 and then running the
regression is in red, and the non-private version in black. The x-axis is the size of the data in log-scale.

5

Sampling from an Inverse-Wishart Distribution (Bayesian Posterior)

In Bayesian statistics, one estimates the 2nd-moment matrix in question by starting with a prior and updating
it based on the examples in the data. More specifically, our dataset A contains n datapoints which we assumed
to be drawn i.i.d from some N (0d , V ). We assume V was sampled from some distribution D over positive
definite matrices, which is the prior for V . We then update our belief over V using the Bayesian formula:
Pr[A | V ]·PrD [V ]
Pr[V | A] = R Pr[A
. Finally, with the posterior belief we give an estimation of V — either by
| W ]·PrD [W ]dW
W
outputting the posterior distribution itself, or by outputting the most-likely V according to the posterior,
or by sampling from this posterior distribution (maybe multiple times). In this work we assume that our
estimator of V is given by sampling from the posterior distribution.
One of the most common priors used for positive definite matrices is the inverse-Wishart distribution.
This is mainly due to the fact that the inverse-Wishart distribution is conjugate prior.6 Specifically, if
6 A family of distributions is called conjugate prior if the prior distribution and the posterior distribution both belong to this
family.

9

our prior belief is that V ∼ Wd−1 (Ψ, k), then after viewing n examples in the dataset A our posterior is
V ∼ Wd−1 (AT A + Ψ), n + k . Here we show that sampling such a positive definite matrix V from our
posterior inverse-Wishart distribution is (, δ)-differentially private, provided the prior distribution’s scale
matrix, Ψ, has a sufficiently large σmin (Ψ). This result is in line with the recent beautiful work of Vadhan
and Zheng [VZ15], who showed that many Bayesian techniques for estimating the means are differentially
private, provided the prior is set correctly. The formal description of our algorithm and its privacy statement
are given below.
Input: A matrix A ∈ Rn×d and a bound B > 0 on the l2 -norm of any row in A.
Privacy parameters:
, δ > 0.
 p

2B 2
Set ψ ←  2 2(n + d) ln(4/δ) + 2 ln(4/δ) .
Sample M ∼ Wd−1 ((AT A + ψ · Id×d ), n + d).
return M and the approximation βe = arg minβd =−1 β T M β.
Algorithm 4: Sampling from an Inverse-Wishart Distribution
Theorem 5.1. Fix  > 0 andδ ∈ (0, 1e ). Fix B > 0. Let
 A be a (n × d)-matrix and fix an integer ν ≥ d.
p
2
2
Let w be such that w = 2B 2 2ν ln(4/δ) + 2 ln(4/δ) /. Then, given that σmin (A) ≥ w, the algorithm
that samples a matrix from Wd−1 (AT A, ν) is (, δ)-differentially private.
We comment on the similarities between Theorem 5.1 and Theorem 3.1. Indeed, the Algorithm 1 essentially samples a matrix from W(AT A + w2 I, k) for some choice of w and k (and then normalizes the
sample by k1 ); and Algorithm 4 samples a matrix from W −1 (AT A + w2 I, k) for a very similar choice of w. In
fact, in Algorithm 1, instead of sampling R and then multiplying it with A, we can sample the same R and
multiply it with
(AT A)1/2 ; or even sample a (r × d)-matrix R̃ where each of its rows is sampled i.i.d from

T
N 0d , A A . (All of those have the same distribution over the output.) And so, much like we did in the
Johnson-Lindenstrauss case, we can also use part of the privacy budget to estimate σmin (AT A) and then set
the parameter ψ accordingly. Details appear in Algorithm 5.
Input: A matrix A ∈ Rn×d and a bound B > 0 on the l2 -norm of any row in A.
Privacy parameters: , δ > 0.
A parameter
 pk0 indicating the minimal
 degrees of freedom.
2
Set ψ ← 4B 2 2k0 ln(8/δ) + 2 ln(8/δ) .
n
o
2
2
Set s ← max 0, σmin (AT A) − 2B ln(2/δ)
+
Z
where Z ∼ Lap( 2B ).

Adjust ψ ← max{0, ψ − s}.
if w > 0 then
Sample M ∼ Wd−1 ((AT A + ψ · Id×d ), k0 ).
else
n
 p

o
2
Set k ∗ ← max k ∈ Z : 4B 2 2k ln(8/δ) + 2 ln(8/δ) ≤ s
Sample M ∼ Wd−1 (AT A, k ∗ ).
end
return M and the approximation βe = arg minβd =−1 β T M β.
Algorithm 5: Sampling from an Inverse-Wishart Distribution whose degrees of freedom are determined
by the input.

10

5.1

Experiments on Sampling from the Inverse Wishart Distribution

To estimate the utility of Algorithms 4 and 5, we conduct similar experiments to before, in the same setting
detailed in Section 3.1. For each choice of  and n we ran 5 predictors. (i) The first one (in black) is the naı̈ve
and non-private Bayesian posterior sampling from the inverse Wishart distrbituion. (ii) The second one is
given by Algorithm 4 (in blue). (iii) The third one is given by Algorithm 5 (in red) where the min-degreesof-freedom parameter is set to n + d (so that we have a direct way to compare between Algorithm 4 and 5).
(iv) The fourth is given by Algorithm 2 where the min-number-of-rows parameter is set 2d (in green), and
(v) the fifth one is Algorithm 5 when the min-degrees-of-freedom parameter is set 2d. This gives us a direct
comparison between Algorithms 2 and 5. We repeat each experiment t = 15 times, measuring the l2 -distance
between the outputted estimator of each algorithm to the true β we used to generate the data. We report the
mean and standard variation of the 15 experiments in Figure 3. The results are again consistent among the
various choices of . Both Algorithm 4 and 5 (techniques (ii) and (iii)) exhibit fairly large errors throughout,
mainly due to the fact that the parameter ψ used in each algorithm depends in n, as opposed to any other
algorithm we present. We were surprised to see how little variance there exists in the results (the variance
is too small to be visible in the figure). We did find it surprising that for the most part, the fact we split
the privacy budget in Algorithm 5 turns out to be consistently costlier than Algorithm 4, even for very large
values of n. Another result that we found interesting is that technique (v) outperforms the JL-technique (iv)
(and it is holds for all values of n). Initially we conjectured that the gap can be explained by Lemma A.1,
where the bound for the inverse-JL has a slightly better second order term than the bound for the standard
JL. However, for some values of n the gap is fairly noticeable, and we leave it as an open problem to see if
this holds for any projection matrix (and not just JL).

6

Comparison to the Analyze Gauss Baseline

In this paper we discuss multiple ways for outputting a differentially private approximation of AT A. One
such way was already given by Dwork et al in their “Analyze Gauss” paper [DTTZ14]. As mentioned
already, Dwork et al simply add to AT A a symmetric matrix N whose entries are sampled i.i.d from a
suitable Gaussian. Furthermore, the magnitude of the noise introduced by the Analyze Gauss algorithm is
the smallest out of all algorithms. Yet, as we stressed before, the output of Analyze Gauss isn’t necessarily a
positive definite matrix. In this work we investigate the effect of these fact on the problem of linear regression.
We study the utility of the Analyze Gauss algorithm for the linear regression problem both theoretically
(the theorem regarding the utility of Analyze Gauss is deferred to Appendix C) and experimentally, in
comparison to the other algorithms we introduce in this work. The high-level message from the experiments
we show here as follows. In the simple case, Analyze Gauss is the best algorithm to use,,7 and when it
returns “unreasonable” answers — so do all other algorithms we use (details to follow). However, there do
exist cases where it under performs in comparison to the additive Wishart noise algorithm (Algorithm 3)
and the Wishart (Algorithm 2) or inverse-Wishart (Algorithm 5) sampling algorithms.
In this section we compare between the following 6 techniques.
1. Analyze Gauss algorithm: output AT A + N with N a symmetric matrix whose entries are i.i.d samples
from a Gaussian (black line, squares.)
2. The JL-based algorithm, Algorithm 2 (blue line, squares.)
3. The additive Wishart noise algorithm given by Algorithm 3 (magenta line, squares.)
4. A scaling version of Analyze Gauss: if the output of Analyze Gauss is not positive definite, add cId×d to
it with c = E[kN k](black line, circles.)
5. Algorithm 5, which, as we commented in the experiments of Section 5, is analogous to Algorithm 2 and
seems to consistently do better than Algorithm 2. Both Algorithm 2 and 5 were given the same min-degreesof-freedom parameter: 2d (blue line, circles.)
6. The scaling version of the additive Wishart random noise, as detailed in the experiment of Section 4.
7 In

our opinion, this result is of interest by itself.

11

Figure 3: (best seen in color) A comparison of the average l2 -error for the estimators based on inverseWishart distribution sampling. The non-private sampler is in black, Algorithm 4 is in blue and Algorithm 5
in red. The JL-based algorithm (Algorithm 2) that effectively samples from the Wishart distribution is in
green; and its analogous algorithm that samples from the inverse-Wishart distribution (Algorithm 5) is in
magenta. The x-axis is the size of the data in log-scale.
√
√
I.e.,
AT A + W − k · V (if this leaves the output positive definite) or AT A + W − ( k − ( d +
p outputting
2 ln(4/δ)))2 · V otherwise (magenta line, circles.)
Post-processing the Analyze Gauss output. We have experimented extensively with multiple ways
to project the output of the Analyze Gauss algorithm onto the manifold of PSD matrices. Indeed, the most
naı̈ve approach is to find a PSD matrix M as to minimize kM − (AT A + N )k2F . Such M effectively turns to
be the result of zeroing out all negative eigenvalue of (AT A + N ). The utility of this approach turns out to
be just as bad as the standard Analyze Gauss algorithm (with no post-processing), returning estimations of
size 12 or 9 when the true β has kβk ≈ 3. Other approached we have experimented with were to try other
values of c for a post-processing of the form AT A + N + cId×d . (Such as setting c to be the upper- and
lower-bound on the singular values of N w.p. ≥ 1 − δ.) The performance of such approaches was, overall,
comparable to the chosen technique (setting c = E[kN k]) but with worse performance then our choice of c.
Therefore, in our experiment, we used the best of all techniques we were able to come up with to post-process
Analyze Gauss. This, however, does not mean that there isn’t another post-process technique for Analyze
Gauss that we didn’t think of which out-performs our own approach.

12

6.1

The Basic Single-Regression Experiment

In the same experiment setting from before (see Section 3.1) we compare our 6 estimators based on the l2 distance to the true β that generates our observations. The results, given in Figure 4 are pretty conclusive:
Analyze Gauss is the better of all algorithms. Indeed, for smaller values of n its output is completely out
of scale (while kβk ≈ 3, the average error of Analyze Gauss is about 9, 12, and sometimes 30). In fact, the
error of Analyze Gauss for small values of n is so large that we don’t even present it in our graphs (and the
standard deviation is so large, that the error bar of Analyze Gauss results in a big spike for such values of n).
However, it is important to notice that for such values of n all other techniques also have a fairly large error
(recall, kβk is roughly 3, so errors > 2.5 essentially give no information about β). Once n reaches a certain
size, then there is a sharp shift transition, and Analyze Gauss becomes the algorithm with the smallest error
for all greater values of n. Eventually, the errors of all algorithms becomes smaller than the error between βb
and β (βb is the non-private estimator of β). We also comment that, like before, technique 5 (Algorithm 5)
in consistently better than technique #2 (Algorithm 2), but also note that both technique have the largest
variances in comparison to all other techniques.

Figure 4: (best seen in color) A comparison of the average l2 -error for our 6 estimators. Analyze Gauss
(squares) and its scaled version (circles) are in black; JL algorithm (squares) and the JL variant that samples
from the inverse Wishart distribution (circles) are in blue; and the additive Wishart noise (squares) and its
scaled version (circles) are in magenta. The x-axis is the size of the data in log-scale.

6.2

The Multiple-Regressions Experiment

In this paper we argue that it is important to use algorithms that inherently output a positive definite
matrix. To that end, we now investigate a more complex case, where the data is close to being singular, such
that additive Gaussian noise is likely to introduce much error. The example we focus on is when the data
A is composed of 2p features: the first p = 20 columns are independent of one another (sampled i.i.d from
13

a normal Gaussian); the latter p = 20 columns are the result of some linear combination of the first p ones.
And so A = [X; y1, . . . , yp ] where for every i we have yi = Xβi + ei where each coordinate of ei is sampled
i.i.d from N 0, σ 2 for σ = 0.5 (fixed for all i). In our experiments, we vary n (from 212 to 227 in powers of
2), but fix  = 0.1. What we also vary is the number of y-features we use in our regression.
Recall, our algorithms approximate the Gram matrix of the data. Once such an approximation is published, it is possible to run as many linear regressions on it as we want — fixing any one column of the data
as a label and any subset of the remaining columns as the features of the problem. This is precisely what
we analyze here. We look at the linear regression problem where the label is some yi0 , and the features of
the problem are the first d columns plus some m additional y-columns.8 (I.e.: {x1 , . . . xp } ∪ {y1 , . . . , ym }
where the latter are disjoint to yi0 .) A good approximation of β should therefore return some βe which is
0 (or roughly 0) on the latter m coordinates. This corresponds to what we believe to be a high-level task
a data-analyst might want to perform: finding out which features are relevant and which are irrelevant for
regression.
The results in this case are far less conclusive and are given in Figure 5. When m = 0, we are back to
the case of a single regression (with no redundant features), and here Analyze Gauss (black, squares) outperforms all other algorithms once n is large enough (in our case, n ≥ 216 ). Yet, it is enough to set m = 1
to get very different results. When m > 0 it is evident that Analyze Gauss really performs badly — in fact,
in most cases its values were far beyond the range of a reasonable approximation for β (taking values like
26 and 45 where kβk ≈ 3.2). The scaled version of Analyze Gauss (black, circles) does perform significantly
better, yet — it is not the best out of all algorithms. In fact, it is consistently worse than the JL-based
algorithms (blue, circles and squares) and from the scaled version of the additive Wishart noise (magenta,
circles) for n < 222 = 4, 194, 304. Note that as m increases, all algorithms’ errors become fairly large. In
addition, Figure 6 shows the variance of our estimators. It is clear that the scaled version of Analyze Gauss
has the smallest variance.9 However, the scaled additive Wishart noise algorithm (magenta, circle) seems to
have a good variance as well, and, as discussed, does out-perform the scaled Analyze Gauss algorithm for a
wide range of values of n.
Discussion. It is possible to interpret the results of this experiment, especially for the larger values of m,
as a detriment for all the algorithms that approximate the Gram matrix of the data. Indeed, we pose the
question of running regression over data where there does exist a large correlation between multiple columns
as an open question. One approach could be to find a differentially private analogues to the techniques
of [Mah11] for choosing a subset of the coordinates that approximate the k-PCA. An alternatively approach
is to analyze the Lasso regression over the output of the algorithms that approximate the 2nd-moment
matrix. In fact, we did experiment (though not extensively) with the Lasso regression. Using off-the-shelf
Lasso regression packages (R package named glmnet), it seems that all algorithms give estimators that are
indeed sparse, but not specifically over the latter m coordinates. Rather, the estimator is sparse both on
the the first p coordinates and on the latter m coordinates. In contract, running the Lasso regression on the
data without additional randomness (non-privately) gives sparsity over the latter m coordinates. We leave
both problems for future work.

References
[BBDS12] J. Blocki, A. Blum, A. Datta, and O. Sheffet. The Johnson-Lindenstrauss transform itself
preserves differential privacy. In FOCS, 2012.
8 We

actually used one more column, of all 1s, representing the intercept.
is a spike up for the largest value of n = 227 , which is recurring throughout our experiments — most estimators
have reasonable performance, but some have a really large error. One possible explanation could be related to the fact that we
shrink sample points to have l2 -norm of at most B.
9 There

14

Figure 5: (best seen in color) A comparison of the average l2 -error for 6 estimators based. Analyze Gauss (squares)
and its scaled version (circles) are in black; JL algorithm (squares) and the JL variant that samples from the inverse
Wishart distribution (circles) are in blue; and the additive Wishart noise (squares) and its scaled version (circles) are
in magenta. The x-axis is the size of the data in log-scale.

[BST14]

Raef Bassily, Adam Smith, and Abhradeep Thakurta. Private empirical risk minimization:
Efficient algorithms and tight error bounds. In 55th IEEE Annual Symposium on Foundations
of Computer Science, FOCS, 2014.

[CMS11]

Kamalika Chaudhuri, Claire Monteleoni, and Anand D. Sarwate. Differentially private empirical
risk minimization. Journal of Machine Learning Research, 12, 2011.

[DFKU13] Paramveer S. Dhillon, Dean P. Foster, Sham M. Kakade, and Lyle H. Ungar. A risk comparison
of ordinary least squares vs ridge regression. JMLR, 14(1), 2013.
[DG03]

Sanjoy Dasgupta and Anupam Gupta. An elementary proof of a theorem of johnson and lindenstrauss. Random Struct. Algorithms, 22(1), January 2003.

[DKM+ 06] Cynthia Dwork, Krishnaram Kenthapadi, Frank McSherry, Ilya Mironov, and Moni Naor. Our
data, ourselves: Privacy via distributed noise generation. In EUROCRYPT, 2006.
[DMNS06] Cynthia Dwork, Frank Mcsherry, Kobbi Nissim, and Adam Smith. Calibrating noise to sensitivity
in private data analysis. In TCC, 2006.
[DRV10]

Cynthia Dwork, Guy N. Rothblum, and Salil P. Vadhan. Boosting and differential privacy. In
FOCS, 2010.

[DS01]

Kenneth R. Davidson and Stanislaw J. Szarek. Local operator theory, random matrices and
banach spaces. In Handbook of the geometry of Banach spaces, volume 1. 2001.

[DTTZ14] Cynthia Dwork, Kunal Talwar, Abhradeep Thakurta, and Li Zhang. Analyze gauss - optimal
bounds for privacy preserving principal component analysis. In STOC, 2014.
[Har13]

Moritz Hardt. Robust subspace iteration and privacy-preserving spectral analysis. In 51st Annual
Allerton Conference on Communication, Control, and Computing, 2013.
15

Figure 6: (best seen in color) A comparison of the standard error of the l2 -error for 6 estimators based,
which are shown in Figure 5. Analyze Gauss (squares) and its scaled version (circles) are in black; JL
algorithm (squares) and the JL variant that samples from the inverse Wishart distribution (circles) are in
blue; and the additive Wishart noise (squares) and its scaled version (circles) are in magenta. The x-axis is
the size of the data in log-scale.

[HK70]

A. E. Hoerl and R. W. Kennard. Ridge regression: Biased estimation for nonorthogonal problems.
Technometrics, 12:55–67, 1970.

[HR12]

Moritz Hardt and Aaron Roth. Beating randomized response on incoherent matrices. In STOC,
2012.

[HTF09]

Trevor J. Hastie, Robert John Tibshirani, and Jerome H. Friedman. The elements of statistical
learning : data mining, inference, and prediction. Springer series in statistics. Springer, 2009.

[KST12]

Daniel Kifer, Adam D. Smith, and Abhradeep Thakurta. Private convex optimization for empirical risk minimization with applications to high-dimensional regression. In COLT, 2012.

[KT13]

Michael Kapralov and Kunal Talwar. On differentially private low rank approximation. In SODA,
2013.

[Mah11]

Michael W. Mahoney. Randomized algorithms for matrices and data. Found. Trends Mach.
Learn., 3(2), February 2011.

[MKB79]

K.V. Mardia, J.T. Kent, and J.M. Bibby. Multivariate analysis. Probability and mathematical
statistics. Academic Press, 1979.

[Sar06]

Tamás Sarlós. Improved approximation algorithms for large matrices via random projections.
In FOCS, 2006.

[Tao12]

T. Tao. Topics in Random Matrix Theory. American Mathematical Soc., 2012.
16

[Tik63]

A. N. Tikhonov. Solution of incorrectly formulated problems and the regularization method.
Soviet Math. Dokl., 4, 1963.

[TS13]

Abhradeep Thakurta and Adam Smith. Differentially private feature selection via stability
arguments, and the robustness of the lasso. In COLT, 2013.

[Ull15]

Jonathan Ullman. Private multiplicative weights beyond linear queries. In Proceedings of the
34th ACM Symposium on Principles of Database Systems, PODS, 2015.

[Upa14]

Jalaj Upadhyay. Differentially private linear algebra in the streaming model.
abs/1409.5414, 2014.

[VZ15]

Salil Vadhan and Joy Zheng. The differential privacy of bayesian inference. Technical report,
Faculty of Arts and Sciences, Harvard University, 2015. Available on http://nrs.harvard.
edu/urn-3:HUL.InstRepos:14398533.

[WFS15]

Yu-Xiang Wang, Stephen E. Fienberg, and Alexander J. Smola. Privacy for free: Posterior
sampling and stochastic gradient monte carlo. In ICML, 2015.

[XKI11]

Bowei Xi, Murat Kantarcioglu, and Ali Inan. Mixture of gaussian models and bayes error under
differential privacy. In CODASPY. ACM, 2011.

A

CoRR,

Useful Lemmas

In this section we detail the main lemmas that we use in our privacy proofs in the following section. The
lemmas and theorems presented here, for the most part, were known prior to our work. We chose to include
so that the uninformed reader can have their full proof, but we make no claim as to the originality of the
proofs of the lemmas. The proofs of Lemma A.1 and Claim A.2 are based in part on the result Dasgupta
and Gupta [DG03] and in part about results regarding the Wishart distribution given in [MKB79] (Theorem
3.4.7). We encourage the reader who is familiar with lemmas and claims in this section to skip their proofs
and turn to Section B where we prove our privacy theorems.
Lemma A.1. Let X be a (r × d)-matrix of i.i.d normal Gaussians (i.e., xi,j ∼ N (0, 1)). Fix β ∈ (0, 1e ).
Then, for any vector v it holds that

 q


2 ln(2/β)
2
Pr v T ( 1r X T X − I)v ≤
2 2 ln(2/β)
+
kvk
≥1−β
r
r
Furthermore, if r ≥ d then denote t =

Pr

T

v (I −

q

2 ln(2/β)
r−d+1

and assume t < 1. Then

1
( r−d+1
X T X)−1 )v


2t − t2
2
≤
kvk ≥ 1 − β
(1 − t)2

(1)

Proof. Fix v. Each entry of Xv is distributed like N (0, kvk2 ) and so v T X T Xv is just the sum of r i.i.d
1
T T
2
Gaussians with variance kvk2 . In other words, kvk
2 v X Xv ∼ χr . Concentration bounds (see Claim A.2)
give therefore that w.p. ≥ 1 − β we have
p
p
√
√
T T
1
( r − 2 ln(2/β))2 ≤ kvk
r + 2 ln(2/β))2
2 v X Xv ≤ (
which implies
 q
−2 2 ln(2/β)
+
r

2 ln(2/β)
r



kvk2 ≤ v T ( 1r X T X − I)v ≤

and so we get the bound on v T ( 1r X T X − I)v.
17

 q
2 2 ln(2/β)
+
r

2 ln(2/β)
r



kvk2

T

v
2
We now argue that v(X vT X)
−1 v ∼ χr−d+1 . To see this, we argue that specifically for the vector ed (the
1
indicator of the d-th coordinate) we have ed (X T X)
∼ χ2r−d+1 , and the results for any v follows from
−1 e
d
taking any unitary function s.t. U T v = kvked , and the observation that the distributions of X and XU T
are identical.
Now, clearly ed (X T X)−1 ed = (X T X)−1
. Now, if we denote the last column of X as xd and the first
d,d



X−d T X−d


d − 1 columns of X as X−d then X T X = 


X−d T xd 
. Thus, the formula for the entries

2
kxd k

xd T X−d

of the inverse give
1
= kxd k2 − xd T X−d (X−d T X−d )−1 X−d T xd
(X T X)−1
d,d

def
= xd I − X−d (X−d T X−d )−1 X−d T xd = xd T P xd
Now, w.p. 1 we have that X−d has full rank (d − 1). For any choice of X−d with full rank we get a matrix P
1
2
which has rank r −(d−1) and its eigenvalues are either 1 or 0. Hence, for any X−d we get (X T X)
−1 ∼ χr−d+1 .
d,d

Since this distribution is independent of X−d we therefore have that this result holds w.p. 1. I.e.:


1
=
z
PDF (X T X)
−1
d,d


Z

T
−1
T
1
= PDF (X T X)
I
−
X
(X
X
)
X
=
P
· PDF I − X−d (X−d T X−d )−1 X−d T = P dP
=
z
−1
−d
−d
−d
−d
d,d

ZP
=


PDFχ2r−d+1 (z) · PDF I − X−d (X−d T X−d )−1 X−d T = P dP

P

Z
= PDFχ2r−d+1 (z) ·


PDF I − X−d (X−d T X−d )−1 X−d T = P dP = PDFχ2r−d+1 (z)

P

Therefore, with probability ≥ 1 − β we have
√

p
p
√
vT v
∈ ( r − d + 1 − 2 ln(2/β))2 , ( r − d + 1 − 2 ln(2/β))2
T
T
−1
v (X X) v
so
√

r−d+1
p
√
r − d + 1 + 2 ln(2/β)

!2
2

kvk ≤ v

T



T
1
r−d+1 X X

−1

√

v≤

r−d+1
p
√
r − d + 1 − 2 ln(2/β)

which implies




1
X T X)−1 v ≤
v T I − ( r−d+1

2

q

2 ln(2/β)
r−d−1

+

2 ln(2/β)
r−d−1

q

ln(2/β) 2
(1 + 2r−d−1
)
q
ln(2/β)
ln(2/β)


2 2r−d−1
− 2r−d−1
1
q
X T X)−1 v ≥ −
v T I − ( r−d+1
ln(2/β) 2
(1 − 2r−d−1
)

Some arithmetic manipulations show that when

2 ln(2/β)
r−d−1

< 1 we have that
q
ln(2/β)
ln(2/β)


2 2r−d−1
− 2r−d−1
T
T
−1
1
q
v I − ( r−d+1 X X)
v ≤
ln(2/β) 2
(1 − 2r−d−1
)
18

!2
kvk2

as this is the larger term of the two.
Claim
√ k and
√ N√(0, 1). Then, for any 0 < ∆ < k we have that
√ let X1 , . . . , Xk be iidPsamples from
P A.2. Fix
Pr[ i Xi2 > ( k + ∆)2 ] < e−∆/2 and Pr[ i Xi2 < ( k − ∆)2 ] < e−∆/2 .
Proof. We start with the following calculation. For any X ∼ N (0, 1) and any s < 1/2 it holds that
sX

E[e

2

Z∞
]=

x2
2
√1 e− 2 esx dx
2π

Z∞
=

x2 (1−2s)
√1 e−
2
dx
2π

√
y=x 1−2s
√
so dy=dx 1−2s

=

−∞

−∞

Z∞

y2
√1 e− 2 √ dy
1−2s
2π

=

√ 1
1−2s

−∞

We now use Markov’s inequality, to deduce that for any λ ∈ (0, 1/2)
2
√
√
√ 2
√
X
Y
P
√
√
2
E[eλ i Xi ]
λXi2 −λ( k+ ∆)2
√
√
Pr[
Xi2 > ( k + ∆)2 ] = Pr[eλ i Xi > eλ( k+ ∆) ] ≤
=
E[e
]e
eλ( k+ ∆)2
i
i

P

k

k
√
√
√
√ 2
2
2
1
2λ
−λ( k+ ∆)2
=
e
= 1+
e−λ( k+ ∆)
1 − 2λ
1 − 2λ


√
√
λk
− λ( k + ∆)2
≤ exp
1 − 2λ


Setting λ =

√
√ ∆√
2( k+ ∆)

Pr[

so that 1 − 2λ =

X

√

√
k√
k+ ∆

we have

 √ √

√
√
√
√
√
Xi2 > ( k + ∆)2 ] ≤ exp λ k( k + ∆) − λ( k + ∆)2

i

= exp

 √
1
2

k∆ −

1
2

√

√
√ 
∆( k + ∆) = exp(− ∆
2)

A similar calculation shows the lower bound.
√
√
√ 2
√ 2
X
Y
P
√
√
2
2
Pr[
Xi2 < ( k − ∆)2 ] = Pr[e−λ i Xi > e−λ( k− ∆) ] ≤
E[e−λXi ]eλ( k− ∆)
i

i

k

√
√ 2
2
1
2λ
=
eλ( k− ∆) = 1 −
1 + 2λ
1 + 2λ


√
√
λk
≤ exp −
+ λ( k − ∆)2
1 + 2λ



Setting λ =

√
√ ∆√
2( k− ∆)



k
2

√

eλ(

√
k− ∆)2

√

so that 1 + 2λ =

√

k√
k− ∆

we have

 √ √

X
√
√
√
√
√
Pr[
Xi2 > ( k + ∆)2 ] ≤ exp −λ k( k − ∆) + λ( k − ∆)2
i

 √
√ √
√ 
= exp − 12 k∆ + 21 ∆( k − ∆) = exp(− ∆
2)

−1
Lemma
√A.3.qFix δ ∈(0, e ). Let X be a matrix sampled from a Wishart distribution Wd (V, m) where
√
m>
d + 2 ln( 2δ ) . Then, w.p. ≥ 1 − δ we have that for every j = 1, 2, . . . , d it holds that



q
√
√
σj (X) ∈ ( m ±
d + 2 ln( 2δ ) )2 σj (V )

19

Furthermore, we also have that for any 0 < α ≤ m it holds


q
√
√
2
kαV − Xk ≤ kV k · |α − ( m −
d + 2 ln( δ ) )2 | and


q
√
√
−1
−1
−1
−1
2
≤ σmin (V ) · |α − ( m +
(αV ) − X
d + 2 ln( δ ) )−2 |

Proof. In order to sample X ∼ Wd (V, m) we first sample a matrix Y ∈ Rm×d in which every entry is i.i.d
normal Gaussian. We then multiply Y by V 1/2 , s.t. every row in Y V 1/2 is sampled i.i.d from N (0d , V ). We
then set X = V 1/2 Y T Y V 1/2 .
Now, we invoke a theorem of Davidson and Szarek [DS01] (Theorem II.13) that states that for any t > 1
we have
√
√
√
√
2
2
Pr[σmax (Y ) > m + d + t] < e−t /2 and Pr[σmin (Y ) < m − d − t] < e−t /2
q
√
√
√
 √
to deduce that w.p. ≥ 1−δ it holds that all of the singular values of Y lie on the interval
m−
d + 2 ln( 2δ ) , m +
d+
Next, we let uj denote the j-th eigenvector of V , corresponding to the j-th eigenvalue σj (V ). Therefore, for
any j we have
T

uj T Xuj = (V 1/2 uj ) Y T Y (V 1/2 uj )




q
q
√
√
√
√
d + 2 ln( 2δ ) )2 kV 1/2 uj k2 = σj (V )( m +
d + 2 ln( 2δ ) )2
≤( m+




q
q
√
√
√
√
uj T Xuj ≥ ( m −
d + 2 ln( 2δ ) )2 kV 1/2 uj k2 = σj (V )( m −
d + 2 ln( 2δ ) )2
and furthermore, for any subspace S we have that
 


q
√
√
max kV 1/2 uj k2
max uT Xu ≤ ( m +
d + 2 ln( 2δ ) )2
u∈S: kuk=1
u∈S: kuk=1

 

q
√
√
min uT Xu ≥ ( m −
d + 2 ln( 2δ ) )2
min kV 1/2 uj k2
u∈S: kuk=1

u∈S: kuk=1

Thus, to complete the first part of the proof, we invoke the Courant-Fischer Min-Max Theorem that
state that
σj (X) =

max

min

{S⊂Rd : dim(S)=j} {u∈S: kuk=1}

uT Xu =

min

max

{S⊂Rd : dim(S)=d−j+1} {u∈S: kuk=1}

uT Xu

Therefore, we can pick S 0 = span{u1 , . . . uj } and S 00 = span{uj , . . . , ud } to deduce


q
√
√
T
2
d + 2 ln( δ ) )2 σj (V )
σj (X) ≥
min
u Xu ≥ ( m −
u∈S 0 :kuk=1


q
√
√
T
2
σj (X) ≤
max
u Xu ≤ ( m +
d + 2 ln( δ ) )2 σj (V )
00
u∈S :kuk=1


As for the second part of the claim, it follows from the fact that αV −X = V 1/2 αI − Y T Y V 1/2 . Now, if we
denote Y = U ΣU T as the SVD decomposition
of Y , q
we have αI
−Y T Y = U (αI − Σ) U T . Since all the entries


√
√
on the diagonal lie in the range |α − ( m ±
d + 2 ln( 2δ ) )2 |. As α ≤ m we have that all eigenvalues are
q


√
√
upper bounded by (m − α) + 2 m
d + 2 ln( 2δ ) and the claim follows. Similarly, for (αV )−1 − X −1 =
q
√


√
V −1/2 αI − Y T Y V −1/2 all eigenvalues lie in the range |α−1 − ( m ±
d + 2 ln( 2δ ) )−2 |, which in this
q
√

√
case is upper bounded by |α−1 − ( m +
d + 2 ln( 2δ ) )−2 |. We comment that the bounds on kαV − Xk
and on k(αV )−1 − X −1 k require we use both the upper- and lower-bounds on the eigenvalues of Y .

20

The other two useful tools we use are the formula for rank-1 updates of the determinant and the inverse
(the Sherman-Morrison lemma).
Theorem A.4. Let A be a (d×d)-invertible matrix and fix any two d-dimensional vectors u, v s.t. v T A−1 u 6=
−1. Then:
det(A + uv T ) = det(A)(1 + v T A−1 u)
(A + uv T )−1 = A−1 −

A−1 uv T A−1
1 + v T A−1 u

Proof. Since we have A + uv T = A(I + A−1 uv T ), we analyze the spectrum of the matrix I + A−1 uv T .
Clearly, for any x ⊥ v we have (I + A−1 uv T )x = x + 0 · A−1 u = x, so d − 1 of the eigenvalues of I + A−1 uv T
1
v, and we have z T (I + A−1 uv T )z =
are exactly 1. As for the last one, take a unit length vector z = kvk
1 + kvk · z T A−1 u = 1 + v T A−1 u. Therefore, det(A + uv T ) = det(A) det(I + A−1 uv T ) = det(A)(1 + v T A−1 u).
As for the Sherman-Morrison formula, we can simply check and see that indeed:
(A + uv T )(A−1 −

B

uv T A−1 uv T A−1
uv T A−1
A−1 uv T A−1
) = I + uv T A−1 −
−
T
−1
T
−1
1+v A u
1+v A u
1 + v T A−1 u


v T A−1 u
1
T −1
−
= I + uv A
1−
=I
1 + v T A−1 u 1 + v T A−1 u

Privacy Theorems

In this section, we provide the formal proofs the our algorithms are differential privacy. We comment that,
because we hope these algorithms will be implemented, we took the time to analyze the exact constants in
our proofs rather than settling for O(·)-notation. In addition to the three algorithms we provide, we give
another theorem about the privacy of an algorithm that adds Gaussian noise to the inverse of the data,
which may be of independent interest.

B.1

Privacy Proof for Algorithm 1

Theorem B.1. Fix  > 0 and δ ∈ (0, 1e ). Fix B > 0. Fix a positive integer r and let w be such that
 q
!

1 + ln(4/δ)
2
2
4
4
2 2r ln( δ ) + 2 ln( δ )
w =B 1+

Let A be a (n × d)-matrix with d < r and where each row of A has bounded L2 -norm of B. Given that
σmin (A) ≥ w, the algorithm that picks a (r × n)-matrix R whose entries are iid samples from a normal
distribution N (0, 1) and publishes R · A is (, δ)-differentially private.
Corollary B.2. assuming  < 1 and δ < e−1 , if it holds that r ≥ 2 ln( 4δ ) then it suffices to have w2 ≥
√
r ln(4/δ)
8B 2
for the results of Theorem B.1 to hold. Alternatively, given input where its least singular value

is publicly known to w, we can set
&
2 '
w2
r=
, if indeed r > 2 ln( 4δ )
8B 2 ln( 4δ )
and satisfy (, δ)-differential privacy. Therefore, if& the rows of A '
are i.i.d draws from a 0-mean multivariate

2
(Σ)
Gaussian with variance Σ, then we may set r as
n σmin
= Ω(n2 ).
4
2
8B ln( δ )

21

Proof. Fix A and A0 be two neighboring (n × d) matrices, s.t. A − A0 is a rank-1 matrix of the form
def

E = A − A0 = ei (v − v 0 )T . We thus denote M as the matrix with the i-th row zeroed out, and have M T M =
AT A − vv T = A0T A0 − v 0 v 0T . Recall that we assume that σmin (A), σmin (A0 ) ≥ w and kEk = kv − v 0 k ≤ 2B.
We transpose A and R and denote X = AT RT and X 0 = (A0 )T RT . For each column yj of RT it holds
that yj T ∼ N (0n , In×n ), and therefore the j-th column of X is distributed like a random variable from
N 0r , AT A . Furthermore, as the columns of R are independently chosen, so are the columns of X are
independent of one another. Therefore, for any r vectors x1 , ..., xr ∈ Rd it holds that
PDFX (x1 , ..., xr ) =
PDFX 0 (x1 , ..., xr ) =

−1
r q
Y

exp − 12 xj T (AT A)−1 xj
(2π)d det(AT A)
j=1
r
Y

q

−1

(2π)d

det(A0T A0 )

exp − 21 xj T (A0T A0 )−1 xj



j=1

We apply the Matrix Determinant Lemma, and the Sherman-Morrison Lemma, and deduce:

det(AT A) = det(M T M ) 1 + v T (M T M )−1 v

det(A0T A0 ) = det(M T M ) 1 + v 0T (M T M )−1 v 0
(M T M )−1 vv T (M T M )−1
1 + v T (M T M )−1 v
(M T M )−1 v 0 v 0T (M T M )−1
= (M T M )−1 −
1 + v 0T (M T M )−1 v 0

(AT A)−1 = (M T M )−1 −
(A0T A0 )−1
Together with the inequality
r
Y
PDFX (x1 , ..., xr )
=
PDFX 0 (x1 , ..., xr ) j=1

s

1+x
1+y

= (1 + x)(1 −

y
1+y )

≤ exp(x −

y
1−y )

for any x, y 6= 1 we have


det(A0T A0 )
exp − 12 xj T ((AT A)−1 − (A0T A0 )−1 )xj
T
det(A A)
1

=


r 
Y
1 + v 0T (M T M )−1 v 0 2
j=1
d
Y

1 + v T (M T M )−1 v

exp − 12 xj T ((AT A)−1 − (A0T A0 )−1 )xj




xj T (M T M )−1 v 0 v 0T (M T M )−1 xj
exp
≤
v (M M ) v −
2
1 + v 0T (M T M )−1 v 0
j=1


1
v T (M T M )−1 v
xj T (M T M )−1 vv T (M T M )−1 xj 
+
−
+
2
1 + v T (M T M )−1 v
1 + v T (M T M )−1 v
P

 

r
xj xj T (M T M )−1 v 0
v 0T (M T M )−1
j=1
1

= exp  r · v 0T (M T M )−1 v 0 −
2
1 + v 0T (M T M )−1 v 0

P
 

r
T
T
−1
T
T
−1
v
(M
M
)
x
x
(M
M
)
v
T
T
−1
j
j
j=1
1
r · v (M M ) v
(2)
· exp  −
+
2
1 + v T (M T M )−1 v
1 + v T (M T M )−1 v
1 

0T

T

−1 0

22

Denote

def

z1 = v 0T (M T M )−1 v 0 − v 0T (M T M )−1  1r

r
X


xj xj T  (M T M )−1 v 0

j=1


def

z2 = v T (M T M )−1 v − v T (M T M )−1  1r

r
X


xj xj T  (M T M )−1 v

j=1

we have that




−z2
r
z1
(v 0 (M T M )−1 v 0 )2
PDFX (x1 , ..., xr )
+
≤
+
ln
PDFX 0 (x1 , ..., xr )
2 1 + v 0 (M T M )−1 v 0
1 + v(M T M )−1 v 1 + v 0 (M T M )−1 v 0

r
≤
|z1 | + |z2 | + (v 0 (M T M )−1 v 0 )2
2
We now turn to analyze each of the above three terms separately. The easiest to bound are the terms
v(M T M )−1 v and v 0 (M T M )−1 v 0 . Weyl’s inequality yields that σmin (M T M ) ≥ σmin (AT A) − B 2 , and so we
 2
−1
2
w
give both terms that bound w2B−B 2 = B
. We turn to bounding |z1 |, |z2 |.
2 − 1
We continue assuming that x1 , . . . , xr were sampled from AT A. If they were sampled from A0T A0 then
the proof is analogous. Denote X as the matrix whose columns are x1 , . . . , xr . We have

z2 = ((M T M )−1 v)T M T M − 1r X T X (M T M )−1 v

= ((M T M )−1 v)T AT A − vv T − 1r X T X (M T M )−1 v



= ((M T M )−1 v)T (AT A)1/2 I − (AT A)−1/2 1r X T X (AT A)−1/2 (AT A)1/2 (M T M )−1 v − (v(M T M )−1 v)2

Recall that X is a matrix whose rows are i.i.d samples from the multivariate Gaussian N 0, AT A .
Therefore, the rows of the matrix X(AT A)−1/2 are i.i.d samples from N (0, Id×d ). In other words, the
distribution of X(AT A)−1/2 is the same as a matrix whose entries are i.i.d samples from N (0, 1). We can
therefore invoke Lemma A.1 and have that w.p. ≥ 1 − δ/2.
q
2
+ 2 ln(4/δ)
) (AT A)1/2 (M T M )−1 v + (v(M T M )−1 v)2
|z2 | ≤ (2 2 ln(4/δ)
r
r
q

≤ (2 2 ln(4/δ)
+ 2 ln(4/δ)
) v T (M T M )−1 (M T M + vv T )(M T M )−1 v + (v(M T M )−1 v)2
r
r

 q

 q
2 ln(4/δ)
2 ln(4/δ)
2 ln(4/δ)
2 ln(4/δ)
T
−1
T
−1 2
+
+ (v(M M ) v) 2
+
+1
= (v(M M ) v) 2
r
r
r
r




 2
−1 q
 2
−2 q
2 ln(4/δ)
w
w
≤ B
2 2 ln(4/δ)
+
+
−
1
2 2 ln(4/δ)
+ 2 ln(4/δ)
+1
2 − 1
2
r
r
B
r
r

23

As the bound on |z1 | is the same as the bound on |z2 | we conclude that



PDFX (x1 , ..., xr )
r
ln
|z1 | + |z2 | + (v 0 (M T M )−1 v 0 )2
≤
PDFX 0 (x1 , ..., xr )
2
−1  p

 2
w
2 2r ln(4/δ) + 2 ln(4/δ)
≤ B
2 − 1

 2
−2  p
3r
w
+ B2 − 1
2 2r ln(4/δ) + 2 ln(4/δ) +
2
!
p
2
2r
ln(4/δ)
+
2
ln(4/δ)
3r

p
+
+ 2
≤

1 + ln(4/δ)
16r ln(4/δ)
(2 2r ln(4/δ) + 2 ln(4/δ))2





1
3
≤
1+
+
<

1 + ln(4/δ)
ln(4/δ) 2 16
by plugging in the value of w2 .

B.2

Privacy Proof for Algorithm 3

Theorem B.3. Fix  ∈ (0, 1) and δ ∈ (0, 1e ). Fix B > 0. Let C1 and C2 be such that they satisfy
√

2 C
√ 2 2 ≤ 2
B
C1 ( C2 − 1)
(E.g., C1 = B 2 and C2 = 14
2 .) Let A be a (n×d)-matrix where each row of A has bounded L2 -norm of B. Let
N be a matrix sampled from the d-dimensional Wishart distribution with ν-degrees of freedom using the scale
matrix V (i.e., N ∼ Wd (V, ν)) for any matrix V with least singular value σmin (V ) ≥ C1 (e.g. V = C1 Id×d )
and ν ≥ bd + 2C2 ln(4/δ)c. Then outputting X = AT A + N is (, δ)-differentially private.
We comment that in order to sample
such an N , one can sample a matrix N 0 ∈ Rν×d of i.i.d normal
√
Gaussians, multiple all entries by B/  and set N 0 = N T N .
Proof. Fix A and A0 that are two neighboring datasets that differ on the i-th row, denoted as v T in A
and v 0T in A0 . Let M denote A or A0 without the i-th row, i.e. M T M = AT A − vv T = A0T A0 − v 0 v 0T .
Therefore, denoting σmin (M ) and σmin (A) as the least singular value of M and A resp., we have that
2
2
2
σmin
(M ) ≤ σmin
(A) ≤ σmin
(M ) + B 2 . Same holds for the least singular value of M and A0 .
Recall that
ν−d−1

PDFWd (V,ν) (N ) ∝ det(N ) 2 exp − 12 tr(V −1 N )
We argue that Wishart-matrix additive noise is (, δ)-differentially private, using the explicit formulation of
the PDF. For the time being, we ignore the issue of outputting a matrix X s.t. either X − AT A, X − A0T A0
or X − M T M are non-invertible. (Note, if our input matrix is A, then Pr[X − AT A non invertible] =
PrN ∼Wd (V,ν) [N non invertible] = 0. However, it is not a-priori clear why we should also have Pr[X −
A0T A0 non invertible] = 0 or Pr[X − M T M non invertible] = 0.) Later, we justify why such events can be
ignored. We now bound the appropriate ratios. If we denote the output of the mechanism as a matrix X,
then we compare

24

PDFWd (V,ν) (X − AT A)
=
PDFWd (V,ν) (X − A0T A0 )

 ν−d−1
2



det(X − AT A)
det(X − A0T A0 )



det(X − M T M − vv T )
det(X − M T M − v 0 v 0T )

=

1

e− 2 (tr(V
 ν−d−1
2

−1

(X−AT A))−tr(V −1 (X−A0 T A0 ))

1

e− 2 (tr(V

−1

(X−AT A−X+A0T A0 T ))

 ν−d−1
2

1 − v T (X − M T M )−1 v
=
exp − 12 tr(V −1 v 0 v 0T ) + 12 tr(V −1 vv T )
0T
T
−1
0
1 − v (X − M M ) v
 ν−d−1

2

1 − v T (X − M T M )−1 v
tr(AB)=tr(BA)
exp − 12 v 0T V −1 v 0 + 12 v T V −1 v
=
1 − v 0T (X − M T M )−1 v 0


We can now use the inequality

ln

PDFAT A+N (X)
PDFA0 T A0 +N (X)

1−x
1−y



= (1 − x)(1 +

y
1−y )

≤ exp(−x +

y
1−y )

for any x and any y 6= 1 to deduce


1 T −1
·v V
− (ν − d − 1)(X − M T M )−1 v
2


1
ν−d−1
T
−1
−1
+ · v 0T
(X
−
M
M
)
−
V
v0
2
1 − v 0T (X − M T M )−1 v 0
≤

Note that we either have X − M T M = X − AT A + vv T = N + vv T or X − M T M = N + v 0 v 0T . And so,
we continue assuming X was sampled using AT A, but the case X was sampled from A0T A0 is symmetric.
Further, we only show a bound for the first term of the two above, as the other term will have the same
upper bound.
T
A)−1 vv T (X−AT A)−1
Note that (X − M T M )−1 = (X − AT A + vv T )−1 = (X − AT A)−1 − (X−A1+v(X−A
, hence
T A)−1 v
(v T (X − AT A)−1 v)2
1 + v T (X − AT A)−1 v
(v 0T (X − AT A)−1 v)2
v 0T (X − M T M )−1 v 0 = v 0T (X − AT A)−1 v 0 −
1 + v T (X − AT A)−1 v
v T (X − M T M )−1 v = v T (X − AT A)−1 v −

=

v T (X − AT A)−1 v
1 + v T (X − AT A)−1 v

And so we have:

v T V −1 − (ν − d − 1)(X − M T M )−1 v

= v T V −1 − (ν − d + 1)(X − M T M )−1 v + 2v T (X − M T M )−1 v

≤ v T V −1 − (ν − d + 1)(X − AT A)−1 v + 2v T (X − AT A)−1 v + (ν − d + 1)(v T (X − AT A)−1 v)2
Now, note that (X − AT A) ∼ Wd (V, ν), and so V −1/2 (X − AT A)V −1/2 ∼ Wd (Id×d , ν). This allows us to
invoke Lemma A.1 to
 −1/2
 !


T
−1/2 −1
V
(X
−
A
A)V
1
v T V −1 − ( ν−d+1
(X − AT A))−1 v = (V −1/2 v)T I −
(V −1/2 v)
ν−d+1

25

and infer that w.p. ≥ 1 − δ/2 we have the following bound

v T V −1 − (ν − d − 1)(X − M T M )−1 v
!
p
2 2(ν−d+1) ln(4/δ) − 2 ln(4/δ)
2
2(ν − d + 1)
√
≤
+ √
+ √
kV −1/2 vk2
√
√
√
( ν−d+1 − 2 ln(4/δ))2
( ν−d+1 − 2 ln(4/δ))2
( ν−d+1 − 2 ln(4/δ))4


p
kV −1/2 vk2
2
2 2(ν − d + 1) ln(4/δ) − 2 ln(4/δ) + 2 +

q
= √
√
( ν−d+1 − 2 ln(4/δ))2
(1 − 2 ln(4/δ) )2
ν−d−1

p
2 2(ν − d + 1) ln(4/δ) − 2 ln(4/δ) + 6 −1/2 2
√
kV
vk
=
√
( ν−d+1 − 2 ln(4/δ))2
Analogously, w.p. ≥ 1 − δ/2 the following bound holds as well:


ν−d−1
T
−1
−1
(X
−
M
M
)
−
V
v0
v 0T
1 − v 0T (X − M T M )−1 v 0

(ν − d − 1)(v 0T (X − M T M )−1 v 0 )2
= v 0T (ν − d − 1)(X − M T M )−1 − V −1 v 0 +
1 − v 0T (X − M T M )−1 v 0

(ν − d − 1)(v 0T (X − M T M )−1 v 0 )2
≤ v 0T (ν − d + 1)(X − M T M )−1 − V −1 v 0 +
1 − v 0T (X − M T M )−1 v 0

(ν − d − 1)(v 0T (X − AT A)−1 v 0 )2
≤ v 0T (ν − d + 1)(X − AT A)−1 − V −1 v 0 +
1 − v 0T (X − M T M )−1 v 0
!
p
2 2(ν−d+1) ln(4/δ) − 2 ln(4/δ)
2(ν − d + 1)
√
≤
+ √
kV −1/2 v 0 k2
√
√
( ν−d+1 − 2 ln(4/δ))2
( ν−d+1 − 2 ln(4/δ))4
Combining the two upper bounds we get




ν−d−1
PDFAT A+N (X)
1 T
−1
T
−1
−
ln
≤ ·v V
(X − M M )
v
PDFA0 T A0 +N (X)
2
1 − v 0T (X − M T M )−1 v 0


ν−d−1
1
(X − M T M )−1 − V −1 v 0
+ · v 0T
2
1 − v 0T (X − M T M )−1 v 0
p
2 2(ν − d + 1) ln(4/δ) − 2 ln(4/δ) + 6 kV −1/2 vk2 + kV −1/2 v 0 k2
√
≤
·
√
( ν−d+1 − 2 ln(4/δ))2
2
p
1
δ< 6
2 2(ν − d + 1) ln(4/δ)
B2
≤
· √
√
σmin (V ) ( ν−d+1 − 2 ln(4/δ))2
All we now need to do is to plug in the fact that ν = bd + C2 · 2 ln(4/δ)c ≥ d − 1 + C2 · 2 ln(4/δ), and
that σmin (V ) ≥ C1 to deduce
√
√


B2
PDFAT A+N (X)
2B 2 C2
2 · 2 ln(4/δ) · C2
p
√
≤
· p
ln
≤
≤
PDFA0T A0 +N (X)
C1 ( C2 · 2 ln(4/δ) − 2 ln(4/δ))2
C1 ( C2 − 1)2

B.3

Privacy Proof for Algorithm 4

Theorem B.4. Fix  > 0 and δ ∈ (0, 1e ). Fix B > 0. Let A be a (n × d)-matrix and fix an integer ν ≥ d.
Let w be such that

 p
B2
w2 =
2 2ν ln(4/δ) + 2 ln(4/δ)

(1 − 2 ln(4/δ) )
26

Then, given that σmin (A) ≥ w, the algorithm that samples a matrix from Wd−1 (AT A, ν) is (, δ)-differentially
private.
We comment on the similarity between the bounds of Theorem B.1 and Theorem B.4. This is after
all quite natural, since the JL-theorem is a way to sample from a Wishart distribution Wd (AT A, r) ( since
every row in the matrix RA is an i.i.d sample from N 0, AT A ). Clearly, one can sample a matrix from
Wd (AT A, r) and invert it, to get a sample from Wd−1 ((AT A)−1 , r) and vice-versa. Therefore, we get similar
bounds. The only slight difference lies in the fact that we require in Theorem B.4 that ν ≥ d, s.t. the
matrix we sample is indeed invertible, whereas we do not require any such lower bound for sampling from
Wd (AT A, r).
Proof. As always, we denote A0 as a neighbor of A that differs just on a single row, which we denote v for
A and v 0 for A0 , and as before, the matrix M is the matrix A with the i-th row all zeroed out. Therefore,
AT A − vv T = A0T A0 − v 0 v 0T = M T M . So, denoting σmin (M ) and σmin (A) as the least singular value of M
2
2
2
and A resp., we have that σmin
(M ) ≤ σmin
(A) ≤ σmin
(M ) + B 2 . Same holds for the least singular value of
0
M and A .
Recall that
ν

PDFW −1 (AT A,ν) (X) ∝ det(AT A) 2 det(X)
d

ν+p+1
2


exp − 21 tr((AT A)X −1 )

y
We invoke the determinant update lemma, the Sherman Morisson lemma and the inequality 1+x
1+y ≤ exp(x − 1+y )
yet again to deduce:

PDFW −1 (AT A,ν) (X)
det(AT A)ν/2 exp − 12 tr((AT A)X −1 )
d

=
PDFW −1 (A0 T A0 ,ν) (X)
det(A0T A0 )ν/2 exp − 12 tr((A0T A0 )X −1 )
d

ν/2


1 + v T (M T M )−1 v
1
T
0T 0
−1
=
exp
−
tr((A
A
−
A
A
)X
1 + v 0T (M T M )−1 v 0
2

 

0T

ν
v (M T M )−1 v 0
1
T
T
−1
T
0 0T
−1
v (M M ) v −
tr((vv − v v )X )
≤ exp
−
2
1 + v 0T (M T M )−1 v 0
2



 1 ν · v 0T (M T M )−1 v 0
1
T
T
−1
T −1
0T −1 0
= exp
ν · v (M M ) v − v X v −
−v X v
2
2 1 + v 0T (M T M )−1 v 0

 


1 0T
ν
1 T
T
−1
−1
T
−1
−1
v ν(M M ) − X
v− v
(M M ) − X
v0
≤ exp
2
2
1 + v 0T (M T M )−1 v 0

We continue assuming X ∼ Wd−1 (AT A, ν) (the case X ∼ Wd−1 (A0T A0 , ν) is symmetric). By definition, we
have that X −1 ∼ Wd ((AT A)−1 , ν). Hence (AT A)1/2 X −1 (AT A)−1/2 ∼ Wd (Id×d , ν), which implies that the
distribution of (AT A)1/2 X −1 (AT A)−1/2 is the same as generating a (ν × d)-matrix of i.i.d N (0, 1) samples
and take its cross-product with itself.
We continue using the Sherman-Morrison formula, and derive the bound


ν · (v T (AT A)−1 v)2
v T ν(M T M )−1 − X −1 v = v T ν(AT A)−1 − X −1 v −
1 − v T (AT A)−1 v


≤ ((AT A)−1/2 v)T νId×d − (AT A)1/2 X −1 (AT A)1/2 ((AT A)−1/2 v)
 p

≤ k(AT A)−1/2 vk2 2 2ν ln(4/δ) + 2 ln(4/δ)

27

which holds w.p. ≥ 1 − δ/2 due to Lemma A.1. Similarly, we have


ν
T
−1
−1
(M
M
)
−
X
v0
− v 0T
1 + v 0T (M T M )−1 v 0

ν · (v 0T (M T M )−1 v 0 )2
= −v 0T ν(M T M )−1 − X −1 v 0 +
1 − v 0T (M T M )−1 v 0

ν · (v 0T (AT A)−1 v)2
ν · (v 0T (M T M )−1 v 0 )2
+
= −v 0T ν(AT A)−1 − X −1 v 0 +
1 − v 0T (M T M )−1 v 0
1 − v 0T (AT A)−1 v 0
0T
T
−1
2

ν · (v (A A) v)
ν · (v 0T (AT A)−1 v)2
≤ −v 0T ν(AT A)−1 − X −1 v 0 +
+
1 − v 0T (M T M )−1 v 0
1 − v 0T (AT A)−1 v 0
 p

≤ k(AT A)−1/2 v 0 k2 2 2ν ln(4/δ) + 2 ln(4/δ)


1
1
+ ν · k(AT A)−1/2 v 0 k2 k(AT A)−1/2 vk2
+
1 − v 0T (M T M )−1 v 0
1 − v 0T (AT A)−1 v 0
Denoting the least singular value of (AT A) as w2 , and using the fact that kvk, kv 0 k ≤ B and crudely upper
bounding v 0T (M T M )−1 v 0 and v 0T (AT A)−1 v 0 by 21 we get
!
 1 B4
PDFW −1 (AT A,ν) (X)
1
B2  p
d
ln
≤ · 2 · 2 2 2ν ln(4/δ) + 2 ln(4/δ) + · 4 (4ν + 4ν)
PDFW −1 (A0T A0 ,ν) (X)
2
w
2 w
d


p
B2
As we have w2 =
2 2ν ln(4/δ) + 2 ln(4/δ) we get that

(1− 2 ln(4/δ) )

ln

PDFW −1 (AT A,ν) (X)
d

PDFW −1 (A0 T A0 ,ν) (X)

!
≤ (1 −


2 ln(4/δ) )

+ 2

4ν
≤
8ν ln(4/δ)

d

B.4

An Additional Privacy Theorem — Gaussian Noise for the Inverse

Theorem B.5. Fix  ∈ (0, 1) and δ ∈ (0, e−1 ). Let A be a (n × d)-matrix where the l2 -norm of each row
T
is bounded by B, where it is publicly known that σminB(A2 A) ≥ 1 + ρ with ρ > 0. Then the algorithm that
outputs (AT A)−1 + N where
N is asymmetric matrix with each entry on or above the main diagonal of N


is sampled i.i.d from N 0, 8 log(2/δ)
ρ2 2

is (, δ)-differentially private.

Proof. The proof of the theorem just bounds the l2 -global sensitivity of the inverse, using the Sherman
Morrison formula. We then use the fact that by independently adding noise to each entry in (AT A)−1
j,k for


2 2 log(2/δ)
is (, δ)-differentially private.
j ≤ k where the noise is sampled i.i.d from N 0, GS2 ·
2
Denote A and A0 , two matrices that differ on a single row, which is denoted v in A and v 0 in A. Therefore,
A A0 = AT A + v 0 v 0T − vv T , so Weyl’s inequality gives that |σmin (AT A) − σmin (A0T A0 )| ≤ B 2 . Denoting
M as the matrix we get by zeroing out the i-th row on A or A0 , we have
0T

(M T M )−1 vv T (M T M )−1
1 + v T (M T M )−1 v
(M T M )−1 v 0 v 0T (M T M )−1
= (M T M )−1 −
1 + v 0T (M T M )−1 v 0

(AT A)−1 = (M T M )−1 −
(A0T A0 )−1
Hence,
0T

0 −1

(A A )

T

−1

− (A A)

T

−1

= (M M )



vv T
v 0 v 0T
−
1 + v T (M T M )−1 v 1 + v 0T (M T M )−1 v 0
28



(M T M )−1

Let x1 , x2 , . . . , xd be the eigenvectors of M , corresponding to the eigenvalues µ1 , . . . , µd . Then, for any j, k
we have



(v 0 · xj )(v 0 · xk )
(v · xj )(v · xk )
−1
xj T (A0T A0 )−1 − (AT A)−1 xk = µ−1
−
µ
j
k
1 + v T (M T M )−1 v 1 + v 0T (M T M )−1 v 0
Due to Weyl’s inequality, σmin (M T M ) ≥ σmin (AT A) − kvk2 ≥ ρ · B 2 . And so, together with the inequality
(a − b)2 ≤ 2a2 + 2b2 we get
(A0T A0 )−1 − (AT A)−1

2
F

=

d
X


(xj T (A0T A0 )−1 − (AT A)−1 xk )2

j,k=1

=

X (v · xj )2 (v · xk )2
2
1 + v T (M T M )−1 v
µ2j µ2k
j,k

X (v 0 · xj )2 (v 0 · xk )2
2
+
1 + v 0T (M T M )−1 v 0
µ2j µ2k
j,k

X
2
≤
(v · xj )2 (v · xk )2 + (v 0 · xj )2 (v 0 · xk )2
2
2
(ρB )
j,k



! 
!
X
X
X
X
X
2
 (v · xj )2 
=
(v · xk )2 +  (v 0 · xj )2 
(v 0 · xk )2 
(ρB 2 )2
j
j
j,k

4

=

C

k

0 4

k

4

4B
4
2kvk + 2kv k
≤
= 2
(ρB 2 )2
(ρB 2 )2
ρ

Utility Theorems

In this section we provide the utility statement for the Analyze Gauss algorithm and the additive Wishart
noise algorithm. Throughout this section we assume our database D ∈ Rn×d is in fact composed of D = [X; y]
where X ∈ Rn×p and y ∈ Rn (so we denote p = d−1). Clearly, to assume y is the last column of D simplifies
the notation, but y can be any single column of D and X can be any subset of the other columns of D.
In this section we will repeatedly use the Woodbury formula, which states that for any invertible A ∈ Rp×p
and U ∈ Rp×k and V ∈ Rk×p of corresponding dimension we have
(A + U V )−1 = A−1 − A−1 U Ik×k − V A−1 U

−1

V A−1

which implies that for any B ∈ Rp×p we have the binomial inverse formula:
(A + B)−1 = A−1 − A−1 (Ip×p − BA−1 )−1 BA−1

(3)

Our goal is to compare the distance between our predictor to the predictor one gets without noise, i.e.
]
T D that approximates D T D, we can decompose it into
to βb = (X T X)−1 X T y. Since we release a matrix D
^
]
^
]
T X and the p-dimensional vector X
T y and compute β
T X)−1 X
T y. We thus give
e = (X
the p × p matrix X
bounds on
^
]
T X)−1 X
T y − (X T X)−1 X T y
βe − βb = (X
Our analysis presents utility analysis that depends on the input parameters. This is in contrast to
previous works on DP ERM that give a uniform bound and obtain it via regularization of the problem.
(This is natural, as for X = 0n×p clearly βb is ill-defined unless we regularize the problem.)

29

Theorem C.1. Fix X ∈ Rn×p and y ∈ Rn s.t. X T X is invertible. Fix η ∈ (0, 1) and ν ∈ (0, 1/e). Denote

^
]
T X = X T X + N and X
T y = X T y + n where each entry of N and n is sampled i.i.d from N 0, σ 2 . Then,
X
√
there exists some constant C ≥ 1 s.t. if we have that σmin (X T X) ≥ 2C
η · σ p log(1/ν), then w.p. ≥ 1 − ν we
have
]
^
T y − (X T X)−1 X T y ≤ 2η β
T X)−1 X
b+ η
βe − βb = (X
C
We comment that this is not precisely the same as the behavior of the “Analyze Gauss” algorithm. The
difference lies in the fact that Analyze Gauss outputs X T X + M where M is a symmetric
matrix whose

entries along and above the main diagonal are sampled i.i.d from a suitable N 0, σ 2 .  However, one can
denote M = √12 (N + N T ) for a matrix N whose entries are i.i.d samples from N 0, σ 2 , and so the same
√
result, up to a factor of 2, holds for Analyze Gauss.
Proof. Plugging in (3) we get

T
−1
^
]
T X)−1 X
Ty = I
(X
(Ip×p − N (X T X)−1 )−1 N (X T X)−1 X T y
p×p − (X X)

+ Ip×p − (X T X)−1 (Ip×p − N (X T X)−1 )−1 N (X T X)−1 n
^
]
T X)−1 X
T y − (X T X)−1 X T y
Denoting Z = (X T X)−1 (Ip×p − N (X T X)−1 )−1 N , we derive a bound on (X
using bounds on kZk, kI − Zk and knk.
√
Standard bounds on a symmetric ensemble of Gaussians [Tao12] give that kN k ≤ C · σ p log(1/ν) w.p.
≥ 1 − ν2 for some suitable constant C > 0. Hence we have that kN k · k(X T X)−1 k ≤ η. Hence, all singular
values of N (X T X)−1 are upper bounded in absolute value by η, and so all singular values of I − N (X T X)−1
η
η
1
lie in the range [1 − η, 1 + η]. This implies that kZk ≤ 1−η
and kI − Zk ≤ 1 + 1−η
= 1−η
. Next we note
p
√
ν
that knk2 ∼ σ 2 · χ2p , and so, w.p. ≥ 1 − 2 it holds that knk ≤ σ( p + 2 ln(2/ν)).
Thus, we get
p
p
2p +
σ
2σ 2 ln(2/ν)
1
η
η
b +
b + η
kβk
·
≤
kβk
βe − βb ≤
1−η
1−η
σmin (X T X)
1−η
C

Corollary C.2. Denote ρ =
we have

^
T X)
σmin (X
√
2σ p log(1/ν) .

Then, for the same constant C in Theorem C.1, if ρ ≥ 2C then
1
2C e
kβk +
βe − βb ≤
ρ
ρ

^
TX
Proof. The proof follows from Theorem C.1, and the observation that we can flip the role of X T X and X
C
because the Gaussian distribution is symmetric. And so, we just use the notation ρ = η .


N n
. Let
Theorem C.3. Let W ∼ Wp+1 (σ 2 I, k), and denote N ∈ Rp×p and n ∈ Rp s.t. W =
nT ∗
X ∈ Rn×p be a matrix s.t. X T X is invertible and let y ∈ Rn , and such that there exists a C ≥ 2 s.t.
p
√
√
^
]
T X = X T X + N and X
T y = X T y + n. Then
σmin (X T X) = C · σ 2 ( k + p + 2 ln(4/ν))2 . Denote X
βe − βb ≤

n p
o
p
√
1
σ 2 (C − 2)
√
b +
kβk
· min 2 2kp · ln(4p/ν), ( k + p + 2 ln(4/ν))2
T
C −1
(C − 1)σmin (X X)

Proof. Because σ 2 I is a diagonal matrix, standard results on the Wishart distribution give that N ∼
Wp (σ 2 Ip×p , k). We therefore denote R as a (k × p)-matrix of i.i.d samples from a normal Gaussian N (0, 1),
and have N = σ 2 RT R. The Woodbury formula gives that
(X T X + N )−1 = (X T X)−1 − σ 2 (X T X)−1 RT (I − σ 2 R(X T X)−1 RT )−1 R(X T X)−1
30

Denoting Q = σR(X T X)−1/2 we get


= (X T X)−1 − (X T X)−1/2 QT (I − QQT )−1 Q (X T X)−1/2
T
T −1
Now, if wedenote Q = U ΛV T where Q’s
 singular values are λ1 , . . . , λd , we get Q (I − QQ ) Q = V ·
λ2i
1
diag 1−λ2 · V T = V · diag 1−λ
· V T . Note that QT Q = σ 2 (X T X)−1/2 RT R(X T X)−1/2 and so,
2 − 1
i
i
i
i
√
√ √
σ 2 ( k+ p+ 2 ln(4/ν))2
due to Lemma A.3 we have λ21 = σmax (QT Q) ≤
≤ C −1 w.p. ≥ 1 − ν/2. Which
σmin (X T X)
1
means that w.p. ≥ 1 − ν/2 we have σmax (QT (I − QQT )−1 Q) ≤ C−1
. And so we have that both (i)
1
C−2
T
−1
T
−1
T
−1
T
−1
(X X) − (X X + N )  C−1 (X X) and (ii) (X X + N )  C−1 (X T X)−1 .
Next we turn to bound knk. One easy bound, given Lemma A.3, is to show that w.p. ≥ 1 − ν/2 it holds
that
p
√
√
knk ≤ kW ed k ≤ kW k · 1 ≤ σ 2 ( k + p + 2 ln(4/ν))2

Alternatively we can derive the following bound. Each coordinate in n is the result of the dot-product
between the j-th column of R, denoted rj with the d-th column of R, denoted rd . Each coordinate in R
is sampled i.i.d from N 0, σ 2 . Next, we use the fact that for two independent Gaussians with the same

)2
)2
variance X, Y ∼ N 0, σ 2 it holds that XY = (X+Y
− (X−Y
with 12 (X + Y ) and 12 (X − Y ) are two
2
2
2
independent10 Gaussians N 0, σ2 . And so rj · rd = Zj1 − Zj2 where Zj1 , Zj2 ∼ √σ2 · χ2k . Tail bounds
for the χ2 -distribution (see Claim A.2) give that w.p. ≥ 1 − ν/2 it holds that each coordinate of n is
p
p
p
2 √
2 √
bounded in absolute value by σ2 ( k + 2 ln(4p/ν))2 − σ2 ( k − 2 ln(4p/ν))2 = 4 2k ln(4p/ν), which
p
means knk ≤ 2σ 2 2k · ln(4p/ν).11
Combining both bounds, we have that w.p. ≥ 1 − ν it holds that

βb − βe = (X T X)−1 − (X T X + N )−1 X T y − (X T X + N )−1 n
2
p
e ≤ 1 k(X T X)−1 X T yk + 2σ (C − 2) k(X T X)−1 k 2kp · ln(4p/ν)
⇒ kβb − βk
C −1
C −1
2
p
1
2σ
(C
−
2)
b +
=
kβk
2kp · ln(4p/ν)
T
C −1
(C − 1)σmin (X X)
p
√
σ 2 (C − 2)
√
e ≤ 1 kβk
b +
or:
kβb − βk
( k + p + 2 ln(4/ν))2
T
C −1
(C − 1)σmin (X X)


 
X +Y
1
is where we need to use the fact that X and Y have the same variance. We have
=
X −Y
1


X +Y
is diagonal iff X and Y have the same variance.
and so the variance of
X −Y
p
11 We conjecture that the true bound in log(p)-factor smaller, i.e. O(σ 2 2kp · ln(4/ν)).
10 This

31

1
−1



X
Y



