Differentially Private Linear Algebra in the Streaming Model
Jalaj Upadhyay
Center for Applied Cryptographic Research
David R. Cheriton School of Computer Science
University of Waterloo.

arXiv:1409.5414v5 [] 25 Oct 2017

jalaj.upadhyay@uwaterloo.ca

Abstract
Numerical linear algebra plays an important role in computer science. In this paper, we initiate the
study of performing linear algebraic tasks while preserving privacy when the data is streamed online.
Our main focus is the space requirement of the privacy-preserving data-structures. We give the first
sketch-based algorithm for differential privacy. We give optimal, up to logarithmic factor, space datastructures that can compute low rank approximation, linear regression, and matrix multiplication, while
preserving differential privacy with better additive error bounds compared to the known results. Notably,
we match the best known space bound in the non-private setting by Kane and Nelson (J. ACM, 61(1):4).
Our mechanism for differentially private low-rank approximation reuses the random Gaussian matrix
in a specific way to provide a single-pass mechanism. We prove that the resulting distribution also
preserve differential privacy. This can be of independent interest. We do not make any assumptions, like
singular value separation or normalized row assumption, as made in the earlier works. The mechanisms
for matrix multiplication and linear regression can be seen as the private analogues of the known nonprivate algorithms. All our mechanisms, in the form presented, can also be computed in the distributed
setting.

Keywords. Differential Privacy, Linear Algebra, Random Projection.

i

1

Introduction

An n × d real-valued matrix is a natural structure for storing data about n entities described by d features.
They arise in many contexts like information retrieval, data analysis, finance, scientific computation, genetics, and machine learning. In most of these applications, one is often required to do various linear algebraic
tasks, like low-rank matrix approximations (LRA), linear regression (L IN -R EG), and matrix multiplication
(M AT-M ULT). For example, LRA is used in latent semantic indexing to speed-up the computation by computing a new representation for each document in the original collection, L IN -R EG is used in finance to
analyze and quantify the systematic risk of an investment, and many matrix operations in scientific computation such as Gaussian elimination, LU decomposition, and the determinant or the inverse of a matrix can
be reduced to M AT-M ULT.
Let us consider some applications of these linear algebraic tasks in learning theory. For example, LRA
is an effective tool in reducing the memory burdens of kernel methods. In kernel methods, we map the data
into a kernel-induced Hilbert space. Given n samples, this requires the calculation of an n × n symmetric,
positive semi-definite kernel matrix, which requires quadratic space and mostly cubic time. This can be very
demanding for large problems, impeding its practical deployment. This is where LRA is helpful. Williams
and Seeger [42] showed that the spectra of kernel matrices decay rapidly. This led to the following approach:
compute the LRA of the kernel matrices and perform computation on those matrices. This simultaneously
brought improvement in many other models, like support vector machine, kernel Fisher discriminant analysis and kernel principal component analysis. Similarly, L IN -R EG, which is used to predict the output for a
new data based on the previous data, has played pivotal role in machine learning. Examples include, capital
asset pricing model which is used to predict demands, supplies, and investment, Reinforcement Learning to
approximate value functions, and learning parameters of a model with linear dynamics (see, for eg., Walsh
et al. [41]). Likewise, M AT-M ULT has found many roles in machine learning other than modelling multivariate regression as a transposed M AT-M ULT problem. For example, Mossel et al. [31] showed how to
improve run time of learning juntas using matrix multiplication.
Motivation of this work. The data in the examples mentioned above and on which we perform the learning contain sensitive information and takes a lot of space. As a result, privacy and memory constraints are as
important as correct computation. The privacy issue has been recently exemplified by the deanonymization
of Netflix datasets, which was cited as one of the reasons to study differentially private low-rank approximation by [22, 23]. Similarly, data used in genetics engineering and finance have natural privacy concerns. This
raises the question of whether one can perform all these tasks while giving a robust guarantee of privacy, like
differential privacy. When privacy is not a concern, there are many algorithms which use almost optimal
space and one-pass over the input matrix (see, for example, Kane and Nelson [27]). On the other hand, all
the known private algorithms that achieve small additive error use multiple passes over the matrix. For example, the traditional Krylov subspace iteration method (on which some of the recent works like [21, 23, 28]
are based) requires a lot of space and multiple passes over the input matrix. The main focus of this work is
to perform linear algebraic computation while preserving privacy of the data with the additional challenge
that the data is received online and we can use sub-linear space.
First sketch based differentially private mechanisms. A natural approach to perform linear algebraic
tasks on online data while preserving privacy is to generate and store a sub-linear space data structure that
can be used to perform these tasks without leaking privacy. In the non-private setting, the standard techniques either perform random sampling or generate a random sketch of the input data. Dwork et al. [14]
showed that known (non-private) random sampling based algorithms for certain statistical queries can be
made differentially private. They also gave an impossibility result for private analogues of sketch based approaches for specific “statistical queries.” This raises doubts over the applicability of sketch based approach
in privacy. In this paper, we show the first set of positive results for sketch based approach on streamed
1

data. We give almost optimal space data-structures that can compute LRA, L IN -R EG, and M AT-M ULT
while preserving differential privacy. Notably, we match the best known space bounds for all these tasks in
the non-private setting. We remark that our results does not contradict the impossibility result of Dwork et
al. [14], which holds for cleverly constructed statistical queries.
P RIVACY MODEL USED IN THIS PAPER . There are two commonly used notions of differential privacy:
event level privacy, where guarantees are at the granularity of individual records in the datasets, and user
level privacy, where guarantees are at the granularity of each user whose data is present in the dataset.
Dwork et al. [13] showed that it is impossible to obtain any non-trivial result with respect to the user level
privacy on a streamed data. Therefore, in this paper, we use event level privacy. More specifically, we say
two matrices A1 and A2 with same rows and columns are neighbouring if the matrix A1 − A2 is a rank-1
matrix and kA1 − A2 k ≤ 1. This notion was also used recently by [22, 23, 21] for differentially private
low-rank approximation and by Blocki et al. [4, 5] and Upadhyay [39] for certain statistical queries. We also
restrict how the mechanism receives the data in the sense that the data matrix is streamed online and arrives
either row-wise or column-wise. We call two data streams neighbouring if they are formed by streaming
entries of neighbouring matrices. With this notion of neighbouring datasets, we use the following definition
for privacy.
Definition 1.1. A randomized algorithm M gives (ε, δ)-differential privacy, if for all neighbouring datastreams A1 , A2 and all S ⊂ Range(M), Pr[M(A1 ) ∈ S] ≤ exp(ε)Pr[M(A2 ) ∈ S] + δ, where the
probability is over the coin tosses of M.
The main focus of this study is the space requirement of the privacy-preserving data-structures.

1.1

P ROBLEM S TATEMENTS AND OUR R ESULTS

In this section, we give the formal description of the problems we investigate in this paper. The performance
of a streaming algorithm is measured by three basic factors: the number of passes over the data stream, the
space used by the data-structure, and the time taken to update the data-structure. All our private mechanisms
for performing linear algebraic tasks are single-pass and achieve almost optimal space bound for one-pass
algorithms. Our private sketch are linear; therefore, our mechanisms extends naturally to turnstile updates.
We reserve the letter n for the number of rows and d for the number of columns of a private matrix. We
assume d < n. For bit complexity, we use the convention of Clarkson and Woodruff [10], i.e., the entries of
a matrix can be represented by κ = log(nd) bits.
L OW- RANK A PPROXIMATION . We start with the problem statement of LRA.
Problem 1. ((α, β, τ )-LRA). Given an n × d matrix A and a target rank k, construct a matrix Ψ with k
orthonormal columns such that kA − ΨΨT AkN ≤ (1 + α) minrank(Ak )≤k kA − Ak kN + τ with probability
at least 1 − β where N represents either Frobenius or spectral norm.
The best known space lower bound in the non-private setting for (α, 5/6, 0)-LRA is Ω(nkα−1 κ) due to
Clarkson and Woodruff [10]. Kane and Nelson [27, Th 6.2] improved the earlier analysis of Clarkson and
Woodruff [10] to show a one-pass non-private streaming algorithm for approximate LRA with respect to the
Frobenius norm with row/column-wise updates where one maintains a data-structure in the form of a sketch
of size O(kα−1 (n + d)κ log(1/β)).
The first approach to prove differential privacy using the sketch generation algorithm of Clarkson and
Woodruff [10] is to use additive noise mechanisms. There are two main drawback of this approach. A
straightforward analysis along the line of [22] with no coherence assumption leads to an additive error, τ , of
order k 3/2 for the Frobenius norm approximation. The second drawback is that it only leads to approximation when the approximation norm is Frobenius. In this paper, we use random Gaussian matrices to generate
the sketch and prove the following in Section 4 (see Theorem 4.1 for the precise statement).

2

Method
Chaudhary et al. [9]

√

Hardt and Roth [22]
Hardt and Roth [23]
Kapralov and Talwar [28]
Hardt and Price [21]
Dwork et al. [17]
This paper (Theorem 4.1(i))
This paper (Theorem 4.1(ii))

Additive noise (τ )
O(nk/ε)
q

kn log(k/δ)
ε

+

µkAkF log(k/δ)
ε

q
2

O( kε

σ1

Privacy Notion
Spectral norm

(rk(A)µ + k log n) log
O(dk 3 /(εγ 2 ))
√


1
δ

log n)

knµ log(1/δ) log(n/γ) log log(n/γ)
εγ 1.5 σk
√
√
O((k n ln(1/δ))/ε) + Õ( k 3 n3/2 /ε2 )

p
O(pnk ln(k/δ)/ε)
O( nk ln(k/δ)/ε)

# Passes
k

Event level

2
√

Event level

k log σ
k
√
k log σ

User level
Event level
Event level

1
1
1

Event level
Spectral norm

Table 1: Comparison Between our Mechanism and Previous Mechanisms for Differentially Private k-Rank
Approximation of an n × d rank-rk(A) matrix A. µ denotes the coherence of A, γ = (σk /σk+1 ) − 1. [22]
and Theorem 4.1(i) bounds the Frobenius norm while the other entries bounds the Spectral norm.

Theorem (Informal). There is an O(kε−1 (n + d)κ) bits data structure that could be used
to publish (ε, δ)-differentially
p private k-rank approximation of an n × dpprivate matrix in a
single pass with τ ≤ O( nk ln(2/δ)/ε) for the Frobenius and τ ≤ O( nk ln(2/δ)/ε) for
the spectral norm.
Note that our data-structure provides LRA in both the spectral and Frobenius norm with the same space
bound as in the non-private setting. The reason behind this is simple. The bound on the single pass algorithm
of [10, 27] and the two-pass algorithm of Sarlos [36] use their bound for M AT-M ULT (see Problem 3), and,
therefore, only achieves a bound on the Frobenius norm. On the other hand, we use perturbation theory
which gives us a unified bound for the both the norms. Notably, our proof does not rely on our bound in
M AT-M ULT (see our result for Problem 3).
C OMPARISON AND TIGHTNESS . Let σ1 , · · · , σrank(A) be the singular values of a matrix A. We compare
our result for LRA with the previous works in Table 1. Due to lack of space, we defer the detail comparison
to Appendix 4.0.1. Here we just state few main points. Works like [9, 17, 28] consider different notion
of neighbouring datasets, so our results are incomparable to theirs. Apart from the result in this paper,
only Dwork et al. [17] give LRA on streamed data though under extra assumptions like normalized row
assumption and a lower bound on the optimal value. Also, we do not make any low-coherence assumption.
Therefore, if we set µ = n (for worst-case) to compare our results, then Table 1 shows that we achieve best
bound on τ with the same privacy model. In fact, we achieve almost optimal τ for spectral norm (see, [23,
Thm. 6.1]). Our result for Frobenius norm is also tight due to the result of Blum et al. [6].
L INEAR R EGRESSION . We start by giving the problem statement.
Problem 2. ((α, β, τ )-L IN -R EG). Given an n × d matrix A and a m set of n × 1 column vectors B =
{b1 , · · · , bm }, output a set of vectors X = {x1 , · · · xm } all in Rd so that for all i ∈ [m], kAxi − bi kF ≤
(1 + α) miny∈Rd×1 kAy − bi kF + τ with probability at least 1 − β/m.
The best known space lower bound for (α, 5/6, 0)-L IN -R EG in the non-private setting is Ω(d2 α−1 κ)
due to [10]. Kane and Nelson [27, Th 6.2] improved the analysis of Clarkson and Woodruff [10] to gave
a (d2 α−1 κ log(1/β)) bits data-structure that can answer (α, β, 0)-L IN -R EG in the non-private setting. We
show that their space bound can be achieved in the private setting at the expense of a small additive error.
Informally, we show the following in Section 5.2 (see Theorem 5.3 for the precise statement).

Theorem (Informal). There is a data-structure that uses O d2 α−1 κ log (1/β) bits of space
√
to compute (α, β, O(s2 α n))-L IN -R EG while providing (ε, δ)-differential privacy, where s =
3

Problem
(α, β, τ )-M AT-M ULT
(α, β, τ )-L IN -R EG
Frobenius norm (α, β, τ )-LRA
Spectral norm (α, β, τ )-LRA

Lower Bound

Ω dα−2 κ

Ω d2 α−1 κ
Ω(nkκ/α)
−

Non-private Algorithm
Õ(dα−2 κ log(1/β))

Õ d2 α−1 κ log(1/β)
Õ(kα−1 (n + d)κ)
−

This paper
Õ(dα−2 κ log(1/β))

Õ d2 α−1 κ log(1/β)
O(kε−1 (n + d)κ)
O(kε−1 (n + d)κ)

Table 2: Our Results with Respect to the Best Known Space Bounds in Non-private Setting.

q

16d log(1/β) ln(2/δ)
αε2

ln



16d log(1/β)
αδ



.

L IN -R EG has been also studied in the local privacy model by Duchi et al. [11] and in the online private
learning model by Jain et al. [25] and Thakurata and Smith [38]. These models are different from ours;
therefore, our results are incomparable to theirs.
M ATRIX M ULTIPLICATION . We start by giving the problem statement.
Problem 3. ((α, β, τ )-M AT-M ULT). An n × d1 matrix A and n × d2 matrix B are given. Output a matrix
C so that kAT B − CkF ≤ αkAkF · kBkF + τ with probability at least 1 − β.
Let d = max{d1 , d2 }. The best known space lower bound for (α, 5/6, 0)-M AT-M ULT in the nonprivate setting is Ω(dα−2 κ) due to Clarkson and Woodruff [10]. Kane and
 Nelson [27, Th 6.2] improved
the analysis of Clarkson and Woodruff [10] to give a O dα−2 κ log (1/β) -bits data-structure that can solve
(α, β, 0)-M AT-M ULT in the non-private setting. In Section 5.1, we show that this bound can be achieved in
the private setting with a small additive error (see Theorem 5.1 for the precise statement).
Theorem (Informal). Given conforming matrices A and B. Let r = O(log(1/β)/α2 ) and s =
p
16r ln(2/δ)ε−1 ln(16r/δ). For large enough d = max{d1 , d2 }, there is data structure which
√
maintains a sketch of size O dα−2 κ log (1/β) , computes (α, β, O(s2 α n))-M AT-M ULT and
provides (ε, δ)-differential privacy.
S PACE OPTIMALITY OF THE DATA STRUCTURES . In Table 2, we give the best known lower bound for
the space required for all the problems we study in this paper, that by the best known data-structures in the
non-private setting, and by our private data-structures. The parameters used in Table 2 are the same as in
Table 1. As one can see, we achieve almost optimal space data-structures for all the three problems. Note
that, we explicitly require Ω to be stored for the mechanism of LRA; in the other two cases, Ω can be picked
from the distribution on the fly. We note that stronger lower bounds on space are achievable for each of
the problems considered because of the non-zero additive error; however, the bounds stated in the table are
presently the best known lower bounds that we can use to make any sort of optimality comparison.
One can also implement our mechanisms as distributed algorithms, a desirable feature as argued by [3].
This is because all the operations used in our mechanisms have efficient distributed algorithms (Jacobi
method for singular value decomposition, Cannon [8]’s algorithm for multiplication, and GMRES [35]’s
residual method).

1.2

O UR T ECHNIQUES

We start by giving an overview of the generic private sketch algorithm and then discuss the ideas used in
this work which are different from other related works.
P RIVATE - SKETCH G ENERATION . In order to get a better bound on the additive noise, we devise a privatesketch generation (PSG) mechanism to generate a private sketch of the input stream. We reduce the privacy
4

guarantee to maintaining a spectral property of the streamed matrix. At a high level, all our mechanisms
use this basic mechanism while maintaining the spectral properties of the private matrix. We note that the
affine transformation of Blocki et al. [4, 5] and Upadhyay [39] could be also used to maintain the required
spectral properties, but this leads to a large additive error. Moreover, as noted by Blocki et al. [4], it does
not guarantee LRA. For these reasons, we use a different method.
M ECHANISM AND ANALYSIS FOR LRA. We follow two-steps to compute a LRA: (i) compute a projection
matrix (range finding step), (ii) compute a LRA by operating the projection matrix on the input matrix
(projection step). In the most naive form, both the steps require the input matrix, and hence two-passes are
required. The first observation is that the mechanism for PSG already gives considerable improvement in the
range finding step. This by itself does not give LRA as the singular values of the projection matrix and private
matrix are not comparable at this stage. We need to use the projection step. The second observation is that,
by a clever use of linear algebra, information gathered in the range-finding step using a random Gaussian
matrix can be used to emulate the projection step without using the input matrix explicitly. However, we
need to reuse the random Gaussian matrix. This complicates the privacy proof (discussed below in more
detail). We also use the idea of Hardt and Price [21] to use an oversampling parameter p to get a sharper
bounds for both the Frobenius and spectral norm.
D IFFERENCES IN TECHNIQUES TO COMPUTE LRA. Our mechanism for LRA is markedly different from
the recent private mechanisms [9, 17, 21, 22, 23, 28]. Among all these mechanisms, only Dwork et al. [17]
computes LRA privately in an online manner by using binary tree technique of Dwork et al. [13]. We differ
at a basic level from these mechanisms: all these works perturb the output by adding noise to it, while we
perturb the input matrix and then multiply noise matrix. Our mechanism uses random Gaussian matrices
(privacy proof holds only for Gaussian matrix); therefore, we use different tricks, like reusing Gaussian
matrix, as mentioned in the last paragraph.
D IFFERENCES FROM PREVIOUS ANALYSES OF LRA. Our analysis has some interesting features. As
mentioned above, we reuse the Gaussian matrix to get a mechanism that uses only one-pass over the input
matrix. In general, reusing randomness can result in privacy breach. Fortunately, the Gaussian matrix is
reused in a specific manner for which we prove that privacy holds under certain spectral property of the
input matrix. We believe this can be of independent interest. Our mechanism works by maintaining the
required singular value for the input matrix. This causes the additive error. The analysis used in this paper
to bound the error differs a lot from the analyses of Clarkson and Woodruff [10] and Sarlos [36]. They use
the trick that a good bound on M AT-M ULT allows a bound on LRA. This limits its applicability to the case
when the approximation metric is the Frobenius norm. Moreover, Sarlos [36] use two-passes over the input
matrix. We use perturbation theory, which allows us to give bounds for both the norms in an unified manner.
We believe this can be of independent interest.
OTHER RESULTS : M ECHANISM FOR M AT-M ULT AND L IN -R EG . At a high level, our mechanisms for
(α, β, τ )-M AT-M ULT and (α, β, τ )-L IN -R EG are private analogues of Clarkson and Woodruff [10]. Their
algorithm use tug-of-war matrices with Rademacher entries (this helps in improving the update time).
Adding Gaussian noise to ensure differential privacy amounts to a large additive error. We use PSG mechanism. For the privacy proof to go through, we need to lift the singular values of the input matrix. Our choice
to lift the singular values is constrained by keeping a check on τ as well as to keep the mechanism one-pass.
In order to use the proofs of Kane and Nelson [27] to give an approximation bound, we use an analogous
variance bound on random Gaussian matrix (since privacy does not necessarily hold for random matrices
used by Kane and Nelson [27]).
R ELATED W ORKS . The first formal definition and mechanism for differential privacy was given by
Dwork et al. [12]. Since then, many mechanisms for preserving differential privacy have been proposed

5

in the literature (see, Dwork and Roth [15]). The literature on non-private streaming algorithms is so extensive that we cannot hope to cover it in any detail here. In the private setting, Dwork et al. [14] studied
pan-privacy, where the internal state is known to the adversary. Subsequently, there have been some works
on online private learning, like Dwork et al. [17], Jain et al. [25], and Thakurata and Smith [38], for various
tasks. There are some recent works on differentially private LRA as well. Blum et al. [6] first studied this
problem and gave a simple “input perturbation” algorithm that adds noise to the covariance matrix. This
was improved by Hardt and Roth [22] who studied LRA for the Frobenius norm under the low coherence
assumption. Kapralov and Talwar [28] and Chaudhary et al. [9] studied the spectral LRA of a matrix by
giving a matching upper and lower bounds for privately computing the top k eigenvectors of a matrix. Hardt
and Roth [23] and Hardt and Price [21] improved their noise bound by proposing robust private subspace
iteration mechanism. Recently, Dwork et al. [17] revisited randomized mechanism to give a tighter bound
and used it to give an private online learning algorithm under a normalized row assumption.

2

Notations and Basic Preliminaries

We reserve the letters A and B for private input matrices, and Ω for a random Gaussian matrix. For an
0
n × d matrix A, we let
 Ai: denote
 the i-th row of A, A:j denote the j-th column of A, and A denote
0 A
the symmetric matrix
corresponding to A. We let At denote the matrix received after t time
AT 0
epochs. The singular value decomposition (SVD) of A is A = VΣUT , where U and V are left and right
eigenvectors of A, and Σ is a diagonal matrix. The entries of Λ are called the singular values of A. Since
U and V are unitary matrices, one can write Ai = VΣi UT for any real value i. We let rank(A) denotes
the rank of the matrix A and σi (A) its singular values. Where it is clear from context, we write σi for the
singular values.
We use various matrix norm. We use the notation k · kFPfor Frobenius norm. A Frobenius norm for a
matrix A = (aij )i∈[n],j∈[d] is defined as following kAkF = ij |aij |2 . For a matrix A, we let kAk2 denote
the 2-norm, i.e., maxx∈Rd kAxk2 /kxk2 . s When we wish to refer to both the Frobenius and the spectral
norm, we overload the symbol k · k and drop the subscript. We let e1 , · · · , ed denote the standard basis
vectors in Rd . We use bold face symbols to denote vectors and 0n to denote an n-dimensional 0-vector. For
a matrix M, we write M  0 if all its eigenvalues are non-zero. In the course of this paper, we use many
standard results from linear algebra and random matrices.
L INEAR A LGEBRA . Our analysis make extensive use of linear algebra and statistical properties of Gaussian distribution. We give an exposition to the level required to understand this paper.
Lemma 2.1. If matrix A and B are conforming, then kABkF ≤ kAk2 kBkF .
Lemma 2.2. Let A and B be Hermittian matrices. If A and B differ in at most one row of Euclidean norm
1, then Tr(AT A) − Tr(BT B) ≤ 2.
U NIVARIATE AND M ULTIVARIATE G AUSSIAN D ISTRIBUTION . A random variable, X, distributed
ac

(x−µ)2
1
cording to a Gaussian distribution has the probability density function, PDFX (x) = √2πσ exp − 2σ2
.
We denote it by X ∼ N (µ, σ 2 ). The Gaussian distribution is invariant under affine transformation, i.e., if
X ∼ N (µx , σx ) and Y ∼ N (µy , σy ), then Z = aX+bY has the distribution Z ∼ N (aµx +bµy , aσx2 +bσy2 ).
The multivariate Gaussian distribution is a generalization of univariate Gaussian distribution. Given
a m dimensional multivariate random variable, X ∼ N
 (µ, Σ), the PDF of a multivariate Gaussian is
1
given by PDFX (x) := √ rank(Σ)
exp − 12 xT Σ† x with mean µ ∈ Rm and covariance matrix Σ =
(2π)

Det(Σ)

E[(X − µ)(X − µ)T ]. If Σ has a non-trivial kernel space, then the PDF is undefined. However, in this paper,
we only need to compare the probability distribution of two random variables which are defined over the

6

On input a streamed column v ∈ R2en of the private matrix, parameters r, ñ, the mechanism samples an r × 2ñ
random matrix Ω, and
PSG1 : Compute Yv = Ωv, and return Yv .
PSG2 : Compute Yv = ΩT Ωv, and return Yv .

Figure 1: Private Sketch Generation (PSG) Algorithm
same subspace. Therefore, wherever required, we restrict our attention to the (sub)space orthogonal to the
kernel space of Σ. Multivariate Gaussian distribution maintains many key properties of univariate Gaussian
distribution. For example, any (non-empty) subset of multivariate Gaussians is a multivariate Gaussian and
linear functions of multivariate Gaussian random variables are multivariate Gaussian random variables, i.e.,
if y = Ax + b, where A ∈ Rn×n is a non-singular matrix and b ∈ Rn , then y ∼ N (Aµ + b, AΣAT ).
We use the following properties of Gaussian matrices, a proof of the first lemma could be found in [37],
while for the rest two, one can refer to [32].
Lemma 2.3. Let Ω be a Gaussian matrix. Then E[kAΩBk2 ] ≤ kAk2 kBkF +kAkF kBk2 for fixed matrices
A and B.
Lemma 2.4. Let Ω be a n × (k + p) Gaussian matrix. Then
s
√
k
e k+p
−1 2
−1
E[kΩ kF ] =
and
E[kΩ k2 ] ≤
.
p−1
p
Lemma 2.5. Let Ω be a random n × (k + p) Gaussian matrix whose entries are picked from the distribution
N (0, 1). Then E[Tr((ΩT Ω)−1 )] = k/(p − 1).
Theorem 2.6. (Johnson-Lindenstrauss lemma) Fix any η < 1/2 and let m be a positive integer. Let Ω be a
k × n matrix with entries picked from a Gaussian distribution N (0, 1), where k ≥ 4(α2 /2 − α3 /3)−1 ln m.
Then for any m unit vector set S in Rn
∀x ∈ S, PrM [kΩxk2 ∈ (1 ± α)kxk2 ] ≥ 1 − 2 exp(−α2 k/8).
D IFFERENTIAL P RIVACY.

We use the following in our analysis explicitly or implicitly.

Theorem 2.7. ([16]). Let ε, δ ∈ (0, 1), and δ 0 > 0. If K1 , · · · , K` are each (ε, δ)-differential private mechanism, then the mechanism K(D) := (K1 (D),
p · · · , K` (D)) releasing the concatenation of each algorithm
is (ε0 , `δ + δ 0 )-differentially private for ε0 < 2` ln(1/δ 0 )ε + 2`ε2 .
Lemma 2.8. Let M (D) be a (ε, δ)-differential private mechanism for a database D , and let h be any
function, then any mechanism M 0 := h(M (D)) is also (ε, δ)-differentially private for the same set of
queries.

3

Differentially Private Sketch Generation

We study differential privacy in the well known streaming model of computation [1]. We present it at the
level required to understand this paper (a more formal definition can be found in [1]). This model has three
entities: a stream generator S, a (database) curator K, and a query maker Q. S starts the process at time
t = 0 and the curator initializes its data structure to D0 . Thereafter, the curator is allowed only one-pass over
the input matrix, i.e., it can access any entry of the data-base during exactly one time epoch, and update its
data structure to Dt using Dt−1 and the newly accessed data-points of the matrix. At any time, t, the query
maker Q makes a query q. The curator responds with q(Dt ). The streaming model has a resource bound on
7

the curator. A curator is only allowed to use total time polynomial in the size of the data-base to construct
the data structure, where the size of data-structure should be sub-linear in the size of the data-base. For
differential privacy, we further require that the response of K to the query of Q should satisfy Definition 1.1
with respect to two neighbouring streams. Differential privacy on streaming data has been also studied as
private online learning by [25] and [38], where the emphasis is on learning about the streamed input – one
bears a regret on a hypothesis evaluated against a data point which is not yet streamed. Due to lack of space,
we do not delve into the detail comparison of the two models. Our main privacy result is as follows.
Theorem3.1. (i) If the singular values
of the matrix whose columns are streamed to PSG1 are at least

p
σmin := 4 r log(2/δ) log(r/δ) /ε, then PSG1 preserves (ε, δ)-differential privacy. (ii) If the singular
values of the matrix whose columns are streamed to PSG2 algorithm are at least σmin := (4r log(r/δ)) /ε,
then PSG2 preserves (ε, δ)-differential privacy.
The proof of part (i) needs some care due to subtleties mentioned later in this section, but the overall
outline follows the idea of Blocki et al. [4]. However, the proof of part (ii) is very different and more
involved. We show that for a streamed matrix A, the PDF of the published matrix using PSG2 is
exp(−Tr(A−1 Φ)/2)∆(Φ)(n−r−1)/2
Q
,
2rn/2 π r(r−1)/4 ∆(A)r/2 ri=1 Γ((n − i + 1)/2)

(1)

P
n
where ∆(·) is the determinant and Φ = ri=1 ai aT
i , where ai ∼ N (0 , In×n ). To achieve equation (2), we
use the chain rule to break the computation to smaller part. We then compute the probability distribution
function of every term in the chain rule. Note that equation (2) is not a Wishart distribution as A−1 Φ need
not be symmetric. This is the most technical part of the proof and we perform this computation from the
basics; rest of the proof requires clever manipulation of the PDF for neighbouring matrices.
PSG2 is interesting in its own regard. It says that two successive applications of a Gaussian matrix in a
defined form preserves differential privacy. In general, reuse of randomness does not preserve privacy, but
what we show here is that if the randomness is reused cleverly, then it is possible to achieve privacy.
Proof. We start with the proof of part (i). We give the proof [4] for the sake of completion. We denote by
e the matrix that differs from A by at most one entry. In other word, if A and A
e differs in row i, then there
A
T
T
T
e
e
e
e
exists a unit vector v such that A − A = E = vei . Let UΣV (UΣV , respectively) be the SVD of A
e respectively).
(A,
e is just a linear transformation of N (0, In×n ).
The PDF for the two distributions, corresponding to A and A,
Therefore,
1
1
PDFAT Y (x) = p
exp(− x(AT A)−1 xT )
d
T
2
(2π) ∆(A A)
1
1 e T e −1 T
PDFA
exp(− x(A
A) x )
e T Y (x) = q
2
e T A)
e
(2π)d ∆(A
We prove the result for a row of the published matrix; the theorem follows from Theorem 2.7. It is
straightforward to see that combination of the following proves differential privacy for a row of published
matrix:
s
h
i
∆(AT A)
e T A)
e −1 ΩxT ≤ ε0 ≥ 1 − δ0 ,
∈ exp(±ε0 ) and Pr xΩT (AT A)−1 ΩxT − xΩT (A
e T A)
e
∆(A
where ε0 = √

ε
4r ln(2/δ)

and δ0 = δ/2r.

8

Q
The first part of the proof follows simply as in [4]. More concretely, we have det(AT A) = i σi2 ,
where σ1 ≥ · · · ≥ σd ≥ σmin are the singular values of A. Let σ
e1 ≥ · · · ≥ σ
ed ≥ σmin be its singular
e Since the singular values of A − A
e and A
e − A are the same, P (σi − σ
value for A.
ei ) ≤ 1 using Linskii’s
i
theorem. Therefore,
!
s
Yσ
X
ei2
ε
p
≤ exp
(e
σi − σi ) ≤ eε0 /2 .
2
σ
32
r
log(2/δ)
log(r/δ)
i
i
i
The second part of the proof is slightly more involved. Each row i of the published matrix is distributed
identically and is constructed by multiplying an n-dimensional vector Ωi that has entries picked from a
normal distribution N (0, 1). Note that E[Ωi: ] = 0n and COV(Ωi: ) = I. Then
h
i
e T A)
e −1 ΩxT = xΩT (AT A)−1 (AT E + ET A)(
e A
e T A)
e −1 ΩxT .
xΩT (AT A)−1 ΩxT − xΩT (A
e =U
eΣ
eU
e T , this simplifies as
Using the singular value decomposition of A = UΣUT and A

 



eΣ
e −1 V
e T )ΩxT .
e T )ΩxT + xΩT (VΣ−2 VT )v eT (U
eΣ
e −2 V
xΩT (VΣ−1 UT )ei vT (V
i
Since x ∼ AT y, where y ∼ N (0, 1), we can further simplify it as


 


eΣ
e −2 V
e T )ΩAT y + yT AΩT (VΣ−2 VT )v eT (U
eΣ
e −1 V
e T )ΩAT y .
yT AΩT (VΣ−1 UT )ei vT (V
i
|
{z
}|
{z
} |
{z
}|
{z
}
t1

t2

t3

t4

e 2 , kΣk2 ≥ w, plugging in the SVD of A and A − A0 = ei vT , and that every term ti in
Now since kΣk
the above expression is a linear combination of a Gaussian, i.e., each term is distributed as per N (0, kti k2 ),
we calculate kti k as below.
k(UΣVT )ΩT (VΣ−1 UT )ei k2 ≤ kΩT k2 ≤ 1,
1
,
σmin
eΣ
e −2 V
e T )Ω(V
eΣ
eU
e T − veT )k2 ≤ kv T (V
eΣ
e −2 V
e T )ΩV
eΣ
eU
e T k2 + kv T (V
eΣ
e −2 V
e T )ΩveT k2
kv T (V
i
i
1
1
≤
+ 2 ,
σmin σmin
eΣ
e −1 V
e T )Ω(V
eΣ
eU
e T − veT )k2 ≤ keT (U
eΣ
e −1 V
e T )Ω(V
eΣ
eU
e T k2 + keT (U
eΣ
e −1 V
e T )ΩveT k2
keT (U
k(UΣVT )ΩT (VΣ−2 VT )vk2 ≤ kΩT k2 kΣ−1 k2 ≤

i

i

≤1+

1

σmin
√

where σmin =

i

i

i

,

r log(2/δ) log(r/δ)
ε



. Using the concentration bound on the Gaussian distribution, each term,


1
t1 , t2 , t3 , and t4 , is less than kti k ln(4/δ0 ) with probability 1−δ0 /2. From the fact that 2 σmin
+ σ21 ln(4/δ0 ) ≤
min
ε0 , we have the second part of the proof.
We now prove that the second variant preserves privacy if the singular values of the streamed matrix
follows the hypothesis of the theorem. We first compute the PDF of the published matrix.
We start by computing the probability density function when the underlying multivariate Gaussian
distribution is N (0, I). The case for arbitrary positive definite covariance matrix follows like the transition from identity to arbitrary positive definite covariance matrices in the multivariate Gaussian distribution. Let g1 , · · · , gr ∼ N (0, I) be r i.i.d multivariate Gaussian distribution, i.e., gij ∼ N (0, 1) for
9

P
1 ≤ i ≤ r, 1 ≤ j ≤ n. The distribution we are interested in is Φ = ri=1 gi giT . We use the notation
PDF(Φ; I) to denote the probability density function of Φ when each random variable is picked using a
normal distribution, i.e., when the covariance matrix of the random variables is I.
Using the chain rule, the joint distribution of the entries of Φ is as follows.
PDF(Φ; I) = PDF(hg1 , g1 i; I)PDF((hg2 , g1 i, hg2 , g2 i)|hg1 , g1 i; I) · · ·
PDF((hgr , g1 i, · · · hgr , gr i)|Φ[r−1] ; I).

(2)

There are r(r + 1)/2 distinct entries, hg1 , g1 i, (hg2 , g1 i, hg2 , g2 i), · · · , (hgr , g1 i, · · · hgr , gr i). Our aim is
to compute each individual term in the product form of the above chain rule. For this, we first analyze and
T = (hg , g i, · · · , hg , g
understand the properties of the distribution of h
that
i 1
i i−1 i). Then we use the fact 
i−1
−1
there is a transformation of Jacobian one from hgi , g1 i, · · · , hgi , g(i−1) i, hgi , gi i − hT
i−1 Φ[i−1] hi−1 to

hgi , g1 i, · · · , hgi , g(i) i to compute each term in the chain rule (see, [32] and [33]).

Distribution of hi−1 . We first prove that hi−1 is an (i − 1)-variate Gaussian distribution. Since the
covariance matrix is I, and g11 , · · · , g1n , · · · , gr1 , · · · , grn are i.i.d. N (0, 1), from the elementary property
of linear functions of normal variables, conditional on gkj for 1 ≤ k ≤ i − 1 and 1 ≤ j ≤ n, hi−1 is
(i − 1)-variate Gaussian distribution with


hg1 , g1 i · · · hg1 , gi i

.. 
..
Φ[i] =  ...
.
. 
hgi , g1 i · · ·

hgi , gi i

Now g11 , · · · , grn , for every j = 1, · · · , n are mutually independent; therefore, we have the covariance


COV(hi−1 , gij ) = COV (hgi , g1 igij ) , · · · , COV hgi , g(i−1) igij = (g1j , · · · , gi−1,j )T
and E[|hi−1 ihhi−1 ||gkj ] = Φ[i−1] for 1 ≤ j < i. This implies
i
h
−1
T
COV hi−1 , gij − hT
Φ
(g
j,
·
·
·
g
)
|g
i−1,j
kj = 0
i−1 [i−1] 1

∀1 ≤ k ≤ i − 1,

(3)

T
as the left hand side equals (g1 j, · · · gi−1,j )T − Φ[i−1] Φ−1
[i−1] (g1 j, · · · gi−1,j ) .

2
P
This implies that hi−1 is independent of kj=1 gij − (g1j , · · · , gi−1,j )T Φ−1
(g
)
. [33]
1j,···
,g
i−1,j
[i−1]
proved that
k 
X

gij − (g1j , · · · , gi−1,j )T Φ−1
[i−1] (g1j , · · · , gi−1,j )

2

∼ χ2n−i+1 ,

(4)

j=1

the standard χ2 distribution with (n − i + 1) degrees of freedom.
Computing every term in the chain rule. From the fact that hi−1 is a (i − 1)-variate Gaussian distribution, equation (3), equation (4), and the identity
∆(Φ[i] ) = ∆(Φ[i−1] )

k 
X

2
gij − (g1j , · · · , gi−1,j )T Φ−1
(g
,
·
·
·
,
g
)
,
i−1,j
[i−1] 1j

j=1

10

where ∆(·) denotes the determinant, we first calculate the joint pdf of


T
−1
hgi , g1 i, · · · , hgi , g(i−1) i, hgi , gi i − hT
i−1 Φ[i−1] hi−1

T

exp − 12 hgi , g1 i, · · · , hgi , g(i−1) i Φ−1
hg
,
g
i,
·
·
·
,
hg
g
i
i
1
i
(i−1)
[i−1]
=
(2π)(i−1)/2 ∆(Φ[i−1 )1/2


(n−i+1)/2−1
−1
hgi ,gi i−hT
i−1 Φ[i−1] hi−1
T Φ−1 h
exp −
hg
,
g
i
−
h
i i
i−1 [i−1] i−1
2
×
2(n−i+1)/2 Γ((n − i + 1)/2)
 T −1

hi−1 Φ[i−1] hi−1
exp −
∆(Φi )(n−i−1)/2
2
= m/2 (i−1)/2
2
π
Γ((n − i + 1)/2)∆(Φ[i−1] )(n−i)/2
= PDF((hgi , g1 i, · · · , hgi , gi−1 i, hgi , gi i)|Φ[i−1] ),

(5)



T Φ−1 h
to
where the last step uses the fact that the transformation from hT
,
hg
,
g
i
−
h
i−1
i
i
i
i−1 [i−1]

hgi , g1 i, · · · , hgi , g(i) i is one-to-one with Jacobian 1.
Computing the joint distribution of Φ. We are done except for putting in the values of every terms in the
chain rule. A simple arithmetic followed by plugging equation (5) in equation (2) gives the closed formed
expression of the pdf of Φ as
!
r
Y
∆(Φ[i] )(n−i−1)/2
exp(−Tr(Φ)/2)∆(Φ)(n−r−1)/2
P
×
Q
2rn/2 π i (i−1)/2 ri=1 Γ((n − i + 1)/2) i=1 ∆(Φ[i−1] )(n−i)/2
=

exp(−Tr(Φ)/2)∆(Φ)(n−r−1)/2
Q
.
2rn/2 π n(n−1)/4 ri=1 Γ((n − i + 1)/2)

Let Φ be the matrix formed in the manner as above with random vectors picked using the distribution
that is defined by N (0n , Σ). Let Φ0 = ΦA. Then we can write the probability distribution function of Φ0
as follows
PDFΦ0 (Φ0 ) ∝ PDFΦ (Φ0 A−1 )

∝ exp Tr −Σ−1 ΦA−1 /2
 
= exp Tr −A−1 Σ−1 Φ /2 ,
where A−1 is the pseudo-inverse of A. This is not distributed like the distribution of Φ.
e be a pair of neighbouring matrices
We can now prove the privacy guarantee. Let δ0 = δ/r. Let A and A
e = E = |vihei | for an unit vector
that form the stream. From our definition of neighbouring matrices, A− A
|vi and some i. The published matrices corresponding to the two neighbouring matrices have the following
probability density function
exp(−Tr(A−1 Φ)/2)∆(Φ)(n−r−1)/2
exp(−Tr(A−1 Φ)/2)
Q
=
C
,
r
2rn/2 π n(n−1)/4 ∆(A)r/2 i=1 Γ((n − i + 1)/2)
∆(A)r/2
e −1 Φ)/2)∆(Φ)(n−r−1)/2
e T A)
e −1 Φ)/2)
exp(−Tr((A
exp(−Tr((A
e I) =
PDF(ΦA;
=
C
,
Q
e r/2 r Γ((n − i + 1)/2)
e r/2
2rn/2 π n(n−1)/4 ∆(A)
∆(A)
i=1
PDF(ΦA; I) =

11


Q
where C = ∆(Φ)(n−r−1)/2 / 2rn/2 π n(n−1)/4 ri=1 Γ((n − i + 1)/2) . It is straightforward to see that
combination of the following proves differential privacy of the published matrix:
s
h 
 
i
e
∆(A)
e −1 Φ | ≤ ε ≥ 1 − δ.
≤ exp(ε/r) and Pr |Tr A−1 − A
(6)
exp(−ε/r) ≤
e
e A)
∆(
Let σ1 ≥ · · · , ≥ σd ≥ σmin be the singular values of A. Let σ
e1 , ≥ · · · , ≥P
σ
ed ≥ σmin be the singular
e
e
e
value for A. Since the singular values of A − A and A − A are the same, i∈G (σi − σ
ei ) ≤ 1 using
Linskii’s theorem, where G is the set of indices for which σi > σ
ei . The first bound follows similarly as in
Blocki et al. [4]. For the second bound required for the privacy, we first bound the following

 

 
e −1 Φ = Tr A−1 (A)(
e A)
e −1 − A
e −1 Φ
Tr A−1 − A

 

 
e −1 − A
e −1 Φ = Tr A−1 EA
e −1 Φ .
= Tr A−1 (A + E)A
e =U
eΣ
eV
e T , and the fact that E = vei T
Using the singular value decomposition of A = UΣVT and A
for some i, we can further solve the above expression.
 
 


e −1 Φ
e −1 Φ = Tr A−1 EA
Tr A−1 − A


eΣ
e −1 V
e TΦ
= Tr VΣ−1 UT vei T U
=

r
X



eΣ
e −1 V
e T gj .
Tr gjT VΣ−1 UT vei T U

j=1

Fix a j. We bound the following.


e e −1 e T
Tr gjT VΣ−1 UT veT
.
i UΣ V gj

(7)

We now look at each term in the above expression. gjT VΣ−1 UT ei as N (0, kVΣ−1 UT ei k2 ), and
eΣ
e −1 V
e T gj as N (0, kvT U
eΣ
e −1 V
e T k2 ). Since v and ei are unit vectors, the norm of the above four
vT U
quantities are less than 1/σmin , and 1 + 1/σmin , respectively.
Therefore, from the concentration inequality of Gaussian distribution, we have




1
1
Pr (7) ≤ 2
+ 2
ln(4/δ0 ) ≤ ε ≥ 1 − δ0 .
σmin σmin

 
e −1 Φ ≤ ε.
Taking union bound, we have with probability 1 − δ, −ε ≤ Tr A−1 − A
Differences between PSG1 and Blocki et al. [4]. At high level, PSG1 has some resemblance to the mechanism of Blocki et al. [4] and Upadhyay [39] if we use random Gaussian matrix as Ω. However, the analogy
ends here, for eg., Blocki
p et al. [4] and Upadhyay [39] perform an affine transformation to convert the private
matrix into a set of { w/n, 1}n vectors, while we perform perturbation to raise the singular values before
invoking PSG (see Sections 4 and 5). The mechanism of [4, Algorithm 3] does not give a guarantee that the
singular values of AT A and their published matrix is close or their eigenvalues are comparable. In other
words, it does not give a LRA. Apart from these major differences, there are couple of subtle differences:
(i) they project the entries of the columns of private matrix to a higher dimensional space; here, we perform
embedding to a lower dimensional subspace in the similar vein as other applications of dimensionality reduction, and (ii) their mechanism uses multiple passes over the input matrix (they require at least two-passes
over the input matrix even with the streaming algorithms for computing the singular value decomposition
(SVD). So, we cannot use their mechanism in any of the problems we study in this paper.
12


On input parameters α, β, ε, δ, the target rank k, set w = ckε−1 ln(k/δ) for a global constant c. Pick a 2n × k
standard Gaussian matrix Ω. On input an n × n matrix A of rank r, the mechanism does the following:


Range Finding. Compute YA = wI A Ω by computing wei Ai: Ω for all streams Ai: , 1 ≤ i ≤ n and
appending them row-wise. Let ΠYA = ΨΨT be the projection matrix corresponding to the range of YA .
Projection: When the whole matrix is streamed, the curator does the following:
1. Let the (unknown) matrix B = ΨT At Ψ. Use the minimal residual method to find a solution to
T
BΨT
t Ω = Ψt Yt .
b t Λt U
b T.
2. Compute the decomposition of Bt = Ūt Λt ŪT , form the product Ût = Ψt Ūt , and publish U
t

t

Figure 2: The Mechanism for k-rank Approximation

4

Low Rank Approximation

Blocki et al. [4] noted that their published matrix is neither close nor their eigenvalues are comparable to
the private matrix. In other words, it does not give a LRA; therefore, we need a different approach. We
use the prototype mentioned in Halko et al. [20], which was also used by Hardt and Roth [22] to improve
the worst-case bound under a low coherence assumption. In this prototype, we construct a low-dimensional
subspace that captures the action of the matrix (range-finding), and then restrict the matrix to that subspace to
compute the required factorization (projection). More concretely, range finding finds a measurement matrix
Y = AΩ, where Ω is a Gaussian matrix in our case and computes the orthonormal projection matrix ΠY
corresponding to the range defined by Y; projection then computes a k-rank matrix B = ΠY A. From this
exposition, it seems that privacy preserving algorithms are required for both the stages; however, we show
that the two-step prototype can be replaced by a two-step algorithm in which the input matrix is explicitly
needed only in the first step at the expense of privacy proof requiring both PSG1 and PSG2 .
We first note that if ΨΨT A is a LRA of A, i.e., kA − ΨΨT Ak ≤ η, then so is ΨΨT AΨΨT . This is
because kA−ΨΨT AΨΨT k = kA−ΨΨT A+ΨΨT A−ΨΨT AΨΨT k ≤ kA−ΨΨT Ak+kΨΨT A−
ΨΨT AΨΨT k ≤ 2η. The crucial observation now for the single pass LRA when Ω is a Gaussian matrix
is that Ω, Y, and the basis for the range of Y contains enough information to compute the matrix B, i.e.,
we do not need A explicitly. The range-finding is private using Theorem 3.1(i), but as we reuse Ω, we have
to rely on Theorem 3.1(ii) to prove the privacy of projection step. In order to simplify the presentation, we
state our mechanism for symmetric matrices in Figure 2 (Hardt and Roth [23] and Kapralov and Talwar [28]
also made this assumption) that computes a LRA of the SVD of A.
The mechanism for LRA presented in Figure 2 assume that the matrix A is provided as the symmetric
rank(A) matrix A0 corresponding to A, and stop updating the data structure once all the rows of A are
streamed. This simplifies the presentation as well as the analysis. By the argument of [22, Fact 2.8], this
leads to a depreciation of the privacy guarantee by half (both ε and δ). The analysis for error bound though
has complications because the right and the left singular vectors of the original matrix might be different.
Keeping this in mind, we present its analysis in the most general form.
Theorem 4.1. Let σ1 ≥ · · · ≥ σrk(A) be the singular values of A. Then for an over-sampling parameter
p with the most common choice being p = k + 1, there is a single-pass mechanism that computes k-rank

13

approximation Ā using O(k(n + d)α−1 κ) bits while preserving (ε, δ)-differential privacy such that
s
1/2

2k
k
(n + d) ln(k/δ)
min kA − A0 kF +
,
and
(i) kA − ĀkF ≤ 1 +
0
p−1
ε
p
rk(A )<k
q
P

1/2
e
(k + p) j>k σj2 2pk(n + d) ln(k/δ)
k
(ii) kA − Āk2 ≤ 1 +
σk+1 +
+
.
p−1
p
ε
The most common choice of p is k + 1, which is what we use in Table 1.
Proof. The space complexity is easy to follow from our convention of bit complexity and because we need
to store the matrix Ω and the sketch. The privacy guarantee follows from Theorem 3.1 and noting that all
the singular-values of the matrix on which Ω is operated from the right is greater than the threshold required
for the statement of the Theorem 3.1, and the distribution of ΨT Ω is the same as that of the second variant
as we reuse Ω (this follows from [7]). Now, it follows from the proof of the second variant (Theorem 3.1)
T from left (or right, respectively) of vectors g , · · · , g ,
that it
not matter if we multiply
1
r
Pr A (or A ) T
Pdoes
r
i.e., i=1 A|gi ihgi | = AΦ and i=1 |gi ihgi |A = ΦAT have the same distribution. Combining all
these arguments, we have the distribution of the second step of projection stage is identical to the second
variant in Figure 1, modulo some deterministic computation. Since, any arbitrary post-processing preserves
differential privacy, we can now complete the proof by invoking Theorem 3.1.
√ The privacy guarantee due
4

k log(2/δ) log(k/δ)

to Theorem 3.1 requires the minimum singular value to be greater than
for the first
ε
k log(k/δ)
variant and
for the second variant. By our choice of w, the singular values of the streamed matrix
ε
√
to the algorithm for PSG are at least the eigenvalues of w2 I + AT A, which are all greater than 16k ln(k/δ)
.
ε
√
4 k log(2/δ) log(k/δ)
Since
 16 k log(k/δ)
, the privacy guarantee follows from Theorem 3.1.
ε
ε
1 and the proof in Theorem 3.1 that, for Φ =
Pr In more detail, it follows from Bura and Pfeiffer [7]
T
i=1 |gi ihgi | and large enough n, the distribution of Ψ Ω is
exp(−Tr(A−1 (D−1 R)(D−1 R)T Φ)/2)∆(Φ(D−1 R)(D−1 R)T )(n−r−1)/2
Q
.
2rn/2 π n(n−1)/4 ∆(A)r/2 ri=1 Γ((n − i + 1)/2)
The proof for the first expression of equation (6) is as before. Following the rest of the steps of Theorem 3.1, using Lemma 2.2, we compute for fixed j ∈ [k + p], instead of equation (7), the second expression
in equation (6) is bounded by the following expression

 

Tr VΣ−1 UT ei vT ṼΛ−2 ṼT αj αjT Tr (D−1 R)(D−1 R)T .
Since R is an orthonormal matrix, using the fact that Tr((DT D)−1 ) = Tr((ΩT Ω)−1 ), Lemma 2.5 with
p = k + 1, and following the remaining steps of Theorem 3.1, we have the privacy result for computing B.
The complete proof follows from the discussion stated at the start of this proof.
Remark 1. We take the liberty to diverge a little to understand the intuition behind Bura and Pfeiffer [7].
One may skip this part without effecting the readability of the rest of this section. The result of Bura and
Pfeiffer [7] uses advance statistical tools, but the intuition can be argued using some basic statistics. On the
other hand, it is well known that the singular values of random matrices are notoriously hard to compute
(see [34]). The basic reasoning behind their proof is the following line of argument. Since the entries of
1

We use the following result of Bura and Pfeiffer [7]: for a random n × k normal matrix Ω, for large enough n, the vector
formed by the entries of its left singular matrix is distributed normally with covariance matrix (D−1 R ⊗ I)(RD−1 ⊗ I)T , where
Ω = LDRT .

14

an n × r matrix Ω is N (0, 1), then for any orthogonal matrices G ∈ Rn×n and R ∈ Rr×r , the entries
of GΩRT is also i.d.d. normal. This can be seen by translating to the vector form of the matrix, i.e.,
v = vec(V) represent rn vector with entries vi+nj = Vij . Then vec(GΩRT ) = (G ⊗ R)vec(Ω). Now
G ⊗ R is also an orthogonal matrix, and multivariate Gaussian distribution is preserved if one multiply by
an orthogonal matrix. Therefore, the distribution of the left singular vectors of Ω is the same as GΩRT .
Consequently, for large enough n, the distribution of each singular vector is also spherically distributed. 2
U TILITY G UARANTEE . In Section 4, we showed that the mechanism for symmetric matrices does what [20,
Section 1.2] prototype algorithm achieves. To construct a mechanism for non-symmetric matrices, we construct two sketches Yt and Ȳt corresponding to A1 and A2 using a single pass over A and usingtwo
Gaussian matrices Ω and Ω̄ of appropriate dimension, where A1 = wI A and A2 = AT wI for
appropriate dimension identity matrices in both A1 and A2 . Basically, we do the following using a single
pass over the matrix A:
  
 
wIn A
Y1
Ω1
Y :=
=
,
Y2
Ω2
AT wId
where In is an n × n identity matrix.
T
Since Y1 Y2 corresponds to a symmetric matrix, we can use the steps used in the projection stage
in Figure 2. Note that Clarkson and Woodruff [10] compute AT A in a single-pass over the matrix A.
Keeping the most general case in mind, we first show that the left singular vectors have hardly any role
T
to play in bounding the perturbation. We assume that we perform SVD.
 Let the SVD of A be UΣV . In
the following
discussion, we compute the approximation of wI A and denote it by A. This is because

wI A is more manageable and any upper bound on the approximation of this matrix is an upper bound
on the approximation of the original matrix. The actual bound as computed in Figure 2 can be obtained
by simply performing the computation on the singular value decomposition as performed in the last step of
mechanism and using the sub-additivity of norms. From the discussion above, we know that Ā = ΠY A;
therefore, we need to bound k(I − ΠY )Ak, where, unless specified, in this section k · k refers to both the
Frobenius as well as the spectral norm. From the Hölder’s inequality on the second moment, we have

1/2
E[k(I − ΠY )AkF ] ≤ E k(I − ΠY )Ak2F
.

(8)

√
We now bound k (I − ΠY )A k. Let Λ = Σ2 + w2 I. Let Λ1 denote the diagonal matrix formed by
the first k singular values and Λ2 be the diagonal matrix for the rest of the singular values. We decompose
VT similarly. Let the matrix formed by the first k rows of VT be V1T and by the rest of the rows be V2T .
Left singular vectors have essentially no role in the approximation bound. Recall that Ω is an 2n × k
T
matrix; therefore, Y = AΩ = U Λ1 V1T Ω Λ2 V2T Ω . It would be useful to consider the first k rows of
Y to be the one that mimics the action of A and the rest of the rows of Y as a small perturbation that we
wish to bound. We first prove that the left singular vectors have essentially no role to play in bounding the
error. Let A0 = UA, then the following chain of equalities are straightforward.
k(I − ΠY )Ak = kUT (I − ΠY )Ak = kUT (I − ΠY )UA0 k = kIA0 − UT ΠY UA0 k.

(9)

Now note that the projection matrix corresponding to a matrix Y is uniquely defined by range(Y), the
range of Y. Therefore, range(UT ΠY U) = UT range(ΠY ) = range(UT ΠY ).
Therefore, kA0 − UT ΠY UA0 k = k(I − ΠUT Y )A0 k. A useful way to understand the above expression
is to view this geometrically and recall that unitary are just rotation in the space: projection by an unitary,
followed by any projection operator, followed by the inverse of unitary brings us to the same space as
projection by an operator followed by the inverse of the unitary.
15

Finding and bounding an appropriate perturbed matrix. We now use the identity that, for two operators O1 and O2 , if the range of O1 is a subset of the range of O2 , then the projection of any matrix
using O1 will have all its norm smaller than the projection by O2 . More concretely, we find a matrix C
such that its range is a strict subset of the range of UT Y. We obtain this matrix by flattening out the first
k rows of UT Y. This is in correspondence with our earlier observation that the first k rows mimic the
action of A and other rows are the perturbation that we wish to bound. Since the first k rows of UT Y
is Λ1 V1T Ω, let C := UT YΩ−1 V1 Λ−1
1 . The rows corresponding to the perturbation are Λ2 V2 Ω. Thus,

−1 T
T
C = I Λ2 V2 V1 Λ1
.
T
Let us denote by S = Λ2 V2T V1 Λ−1
1 . It is not difficult to see that range(C) ⊂ range(U Y). Moreover,
ΠC  I, ΠUT Y ΠC ΠUT Y  ΠUT Y . This follows from the fact that range(C) ⊂ range(UT Y) and the
following derivation
ΠUT Y  ΠUT Y ΠC ΠUT Y = ΠC ΠUT Y = (ΠUT Y ΠC )T = ΠC .
An immediate result of the above is the following:
k(I − ΠUT Y )A0 k ≤ k(I − ΠC )A0 k.

(10)

Since, ΠC = C(CT C)−1 CT , we have the following set of derivations.
 
 
 
−1


 I −1
I 
I
T
T
I + ST S
I ST
I S =
I S
ΠC =
S
S
S




T
−1

(I + ST S)−1
(I + ST S)−1 ST
I(I + S S)
T
I S =
=
S(I + ST S)−1 S(I + ST S)−1 ST
S(I + ST S)−1


(I − ST S)
(I + ST S)−1 ST
,

S(I + ST S)−1
0
where the last inequality uses the fact that I − ST S  (I + ST S)−1 and S(I + ST S)−1 ST  0. Therefore,


ST S
I − (I + ST S)−1 ST
.
I − ΠC 
I − ((I + ST S)−1 ST )T
I

Conjugating I−ΠC with Λ, and applying the fact that for every positive definite matrix, P =


X Y
,
YT Z

we have kP k ≤ kXk + kZk, we get
k(I − ΠC )A0 k ≤ kST SA0 k + kA0 k
for any norm. From here on, it is easy arithmetic to show that
q
0
k(I − ΠC )A k ≤ kΛ02 k + kΛ02 V2T Ω(V1T Ω)−1 k

(11)

(12)

for both the required norms. Using equation (9), equation (10) and equation (11), this gives us a bound on the
approximation of matrix A0 . Till this point, our analysis closely follows the ideas of [20], accommodating
the steps of our algorithm. Now, all that remains is to bound kΛ2 V2T Ω(V1T Ω)−1 k, and for this, we have to
analyze the matrix Ω.

16

E RROR BOUND FOR F ROBENIUS NORM We now exploit the rotational invariance of a Gaussian distribution. An important point to note is that (ΩΩT )−1 exists and has a well defined trace. The first part of
−1
the right hand side of equation (12) is immediate. Thus, if we bound E[kΛ02 VT
2 Ω(V1 Ω) k], we are done.
This could be accomplished as below.
v 

u
u X
u 
−1
E[kΛ2 VT
|(Λ2 )0 Πij (V1 Ω−1 )jj |
2 Ω(V1 Ω) k] ≤ tE
ij

ij

q
≤ kΛ02 kF kΩ−1 kF
r

 q
= kΛ02 kF Tr (ΩΩ−1 )T = kΛ02 kF Tr(ΩΩT )−1
q
q
≤ Tr(ΩΩT )−1 min kA − A0 kF + (n + d)wTr(ΩΩT )−1 .
rk(A0 )≤k

The utility guarantee follows by plugging this value in equation (12), and combining equation (8) and
the fact that (ΩΩT )−1 has a well defined trace k/(p − 1) [32].
E RROR BOUND FOR S PECTRAL NORM In order to bound the second term, we use few well known
facts in the theory of random matrices to simplify equation (12). In particular, using Lemma 2.3 and 2.4,
and Holder’s inequality, the statement of the theorem for the spectral norm follows. The utility bound then
follows using the same arithmetic of representing Λ0 in terms of Λ as done in the case of Frobenius norm.
In more details, we first bound
h
i
 h
i1/2
 h
i
E kΛ02 V2T Ω(V1T Ω)−1 k ≤ kΛ02 k E k(V1T Ω)−1 k2F k
+ kΛ02 kF E k(V1T Ω)−1 k

1/2


≤ kΛ02 k E kΩ−1 k2F k
+ kΛ02 kF E kΩ−1 k ,
and then invoke Lemma 2.4 followed by the sub-additivity of norms. Making these substitution and on
simplification, we get the bound stated in Theorem 4.1(ii).
q
P
p

1/2
e
(k + p) j>k λ2j
2 k(n + d) ln(2/δ)
k
kΛ2 k2 +
+
.
(12) ≤ 1 +
p−1
p
ε
4.0.1

Tightness of Bounds and Comparison with Earlier Works

We compare our results with the best possible results in the non-private setting. Eckart and Young [18] have
shown that the quantity minrk(A0 )<k kA − A0 kF in the first term of Theorem 4.1(i) is optimal. Likewise,
Mirsky [29] proved that λk+1 is the minimum spectral error for k-rank approximation. The second term
in Theorem 4.1(ii) shows that we also pay for the Frobenius norm error when doing a unified analysis.
However, when the oversampling parameter p ≈ k, then the factor on λk+1 is constant and that on the
second term is of order k −1/2 . In fact, on closer analysis,
!
p
p
e( (k + p) min{d, n} − k)
2 k(n + d) ln(2/δ)
k
+
λk+1 +
,
kA − Āk2 ≤ 1 +
p−1
p
ε
therefore, the error always lies within some polynomial factor of λk+1 , modulo some additive error. As
pointed out of Halko et al. [20], one can improve this by power-iteration, the method used by Hardt and
Price [21] and Hardt and Roth [23] and multiple pass mechanism. However, it seems unlikely to improve it
in a single-pass.
17

Kapralov and Talwar [28] showed a lower bound on additive error when δ = 0 for neighbouring data
differing by unit spectral norm. Our privacy proof depends strongly on the fact that δ 6= 0. In fact, our
bound is vacuous if δ = 0. Though incomparable due to difference in the notion of neighbouring data, this
separation gap further strengthen the belief that better bounds are possible for δ 6= 0. Recently, Dwork et
al. [17] also showed a bound in the online learning model, which is a different model of computation.
We compare our results with the works stated in Table 1 in more detail.
Chaudhary et al. [9]: They give low-rank approximation in the spectral norm. Additionally, they achieve
(ε, 0)-differential privacy, which is only achieved by Kapralov and Talwar [28]. Their definition of
neighbouring data sets can be (arguably) considered the most general in the sense that they consider
two data sets neighbouring if they differ by at most one in the spectral norm. They use exponential distribution to sample a singular vector and give a heuristic, but practical implementation using Markov
chain Monte-Carlo. On the negative side, their mechanism uses k rounds; therefore, it cannot be implemented in a streaming fashion. Since the notion of neighbouring data-sets and privacy guarantee
achieved is different from that of ours, we believe our result is incomparable to that of Chaudhary et
al. [9]. However, if we just concentrate√on the additive error bound, they achieve a bound of order
O(nk/ε) compared to our bound τ ≤

2

k(n+d) ln(2/δ)
ε

(Theorem 4.1(ii)).

Hardt and Roth [22]: The authors use two passes over the input matrix; therefore, it does not fall in our
one-pass streaming model of computation. They use the same notion of neighboring data-sets as we
do in this paper. This makes their coherence conditions and notion of neighbouring data sets rotationally invariant. As argued by Blocki et al. [4], we achieve a better utility bound in the range
finding step. Intuitively, this could be seen as a consequence of the absence√of additive Gaussian noise. More concretely, Hardt and Roth [22] achieved an error bound of kn log(k/δ)/ε
p
√+
µkAkF (n/d)1/2 log(k/δ)/ε. Their error bound depends on kAkF , which can be as large as nd
for binary matrices in the worst case when the matrix is not as well-behaved as captured by lowcoherence assumption. On the other hand, we achieve a bound that is independent of kAkF .
Hardt and Roth [23]: In some sense, this paper is based on Krylov subspace iteration combined with powering method of Halko et al. [20]. They define two data-sets as neighbouring in the same manner
as in Hardt and Roth [23]. The coherence definition used in this paper depends on the maximum
value of the left or right singular vectors, and is, therefore, rotationally variant. This work assumes
that the singular value are well separated, i.e., the first and the k-th singular value has a non-trivial
separation, and give LRA in spectral norm. Their bound, however, depends on the rank of the input matrix. Their mechanism uses k rounds of subspace generation, each of which depends on the
spectrum of the matrix and uses the power-iteration method [20]. This helps them in achieving better multiplicative bound, but make them unsuitable in a streaming model. A note on multiplicative bound is due here. We believe that the general application of LRA is for thin matrices with
very small tail singular values (for example, see, [42]). Therefore, we feel that in practical scenario,
polynomial multiplicative error would not be that big an issue. The additive error bound
√ computed
p
2
k(n+d) ln(2/δ)
by [23] τ ≤ O(k 2 ε−1 (rank(A)µ + k log n) log(1/δ) log n) compared to τ ≤
ε
(Theorem 4.1(ii)).
Kapralov and Talwar [28]: The only assumption this paper makes is that of singular value separation of
the same form as in Hardt and Roth [23]. They also give low-rank approximation in the spectral norm.
Additionally, they achieve (ε, 0)-differential privacy, which is only achieved by Chaudhary et al. [9].
They use the same definition as in Chaudhary et al. [9]. They also sample a singular vector from exponential distribution, but they give a net-based algorithm to perform the sampling in polynomial time.

18

On the negative side, their mechanism uses k rounds; therefore, it cannot be implemented in a streaming fashion. Since the notion of neighbouring data-sets and privacy guarantee achieved is different
from that of ours, we believe our result is incomparable to that of [28]. However, if we just concentrate on the additive error bound, they achieve a bound of √
O(dk 3 /(εγ 2 )), where γ is the separation
between the singular values, compared to our bound τ ≤

2

k(n+d) ln(2/δ)
ε

(Theorem 4.1(ii)).

Hardt and Price [21]: In this recent work, Hardt and Price [21] gave a robust subspace iteration mechanism that allows to publish LRA with noise independent of the rank of the input matrix, thereby,
resolving one of the open problems in [23]. They define neighbouring data-sets in the same manner
as in [22, 23]. However, they also make an assumption on the singular value separation–a separation between the k-th and (k + 1)-th singular value of the input matrix. Their mechanism uses k
rounds of subspace generation, each of which depends on the spectrum of the matrix to reduce the
multiplicative
error; therefore, it cannot be implemented in a streaming fashion. We achieve a bound
√
p
2 k(n+d) ln(2/δ)
of
(Theorem 4.1(ii)) compared to λ1 knµ log(1/δ) log(n/γ) log log(n/γ)/εγ 1.5 λk
ε
of [21].
Dwork et al. [17]: Dwork et al. [17] gave the first single-pass online learning algorithm for private lowrank approximation under the assumption that the rows of the input matrix are normalized. They
consider the online-learning model, which is very different from our model, and we do not see any
natural way to compare. They use the follow the perturbed leader (FTL) algorithm of Kannan and
Vempala [26] with the binary tree technique of Dwork et al. [13]. This idea was previously used by
√
Jain et al. [25] as well. They give a bound that assumes a lower bound of k n log2 (m/δ)/ε2 on
the optimal value,√where δ < 1/ poly(n). More concretely, if OPT is the optimal value, then their
error bound is O( kOPTn1/4 log2 (m/δ)). We do not make any of the assumptions made by them
√
and, if we just consider the end result, we achieve a bound which is factor k n better than theirs
(see Theorem 4.1(ii)). The case that we are able to bypass their lower bound gives a mathematical
indication that the unit norm notion of neighbouring data is strictly weaker than user-level privacy.

5

Other Applications of Private Sketch Generation

In this section, we give two applications of PSG1 : mechanisms for M AT-M ULT and L IN -R EG. Theorem 3.1
guarantees privacy if the singular values are high enough, so our basic approach would be to lift the singular
values of the private matrix above the threshold of Theorem 3.1(i). However, we have to be careful. For
example, if we use the affine transformation based approach of Blocki et al. [4] and Upadhyay [39], then it
would lead to an additive error proportional to the Frobenius norm of input matrices. To control the additive
error, we follow a different approach. We first transform any conforming matrices A and B to A0 and B0 ,

T

respectively, and then use the identity, I A0 I B0 = I + A0 B0T to perturb the input matrix with a
careful choice of parameters. Intuitively, τ is due to the identity term of the published matrix. We use the
same idea for L IN -R EG as well.

5.1

Matrix Multiplication

We present the mechanism for M AT-M ULT below and the result is stated in Theorem 5.1. The main idea is
to lift the spectra of the input matrices above the threshold of Theorem 3.1(i).
I NITIALIZATION . On input parameters α, β, ε, δ, set r = O(log(1/β)/α2 ). Set s =
intial sketches of A and B to be all zero matrices YA0 and YB0 .

19

q
16r ln( 2δ )/ε ln( 16r
δ ). Set the

DATA - STRUCTURE UPDATE . Set d = max{d1 , d2 }. On input a column a of an n × d1 matrix A and column

b :a = sea 0n+d A:a and B
b :b =
b of an n × d2 matrix B at time epoch t, set the column vector A

n+d
b
b
seb 0
B:b . Invoke PSG1 with inputs (A:a , r, n+d) and (B:b , r, n+d). Update the sketches by replacing the columns a of YAt−1 and b of YBt−1 by the respective returned sketches to get the sketch YAt , YBt .
T
A NSWERING MATRIX PRODUCT. On request to compute the product at time t, compute YA
YBt .
t

Theorem 5.1. Let Ω be the random matrix used by PSG. Then, the data-structure generated by mechanism above uses
O(dα−2
√κ log(1/β)) bits of space, and on input conforming matrices A and B, computes (α, β, τ )-M AT-M ULT with
τ ≤ s2 nα additive error and (ε, δ)-differential privacy.
Proof. The proof of the utility of mechanism in Section 5.1 follows readily from Lemma 5.2, which is the variance
bound computed by [27].
Lemma 5.2. Let Ω be a r × n matrix as constructed in Section 3 with every entries picked from the distribution
N (0, 1), then for a set of m vectors, v1 , · · · , vm ∈ Rn , with probability at least 1 − 2 exp(−rα2 /8), for any pair
vi , vj , we have |hΩvi , Ωvj i − hvi , vj i| ≤ αkvi k · kvj k.
The basic intuition is that the multiplication of the scaled identity matrix causes the additive error while the
multiplicative error is due to the result of [27]. We follow up with the details. We need to upper bound the quantity
b T ΩT ΩB/s
b 2 kF . First note that
kAT B − A


b T ΩT ΩB
b = sI 0 0 AT ΩT Ω sI 0 0 BT T = (s2 IΩT ΩI + AT ΩT ΩB).
A
Therefore,
b T ΩT ΩBk
b F = kAT B − AT ΩT ΩB − s2 IΩT ΩIkF
kAT B − A
≤ kAT B − AT ΩT ΩB + s2 I − s2 IΩT ΩIkF
≤ kAT B − AT ΩT ΩBkF + s2 kI − ΩT ΩkF .

(13)

To bound the first term, let random variable Xij denote (AT B)ij − (AT ΩT ΩB)ij . Then, with probability at least
1 − 2 exp(−kα2 /8), we have |Xij | ≤ αkA:i k2 · kB:j k2 . Using Lemma 5.2, this results in
X
X
2
kAT B − AT ΩT ΩBkF =
|Xij |2 ≤
α2 kAi: k22 kB:j k22 ≤ α2 kAk2F kBk2F .
(14)
For the second term, we need to bound the variance on unitaries. This follows from the following set of inequalities.
T
T
T
kUT
1 Ω ΩU2 − U1 U2 k2 = kU1 (Ω Ω − I)U2 k2 = kΩ Ω − Ik2


xT (ΩT Ω − I)x
= max
≤ ((1 + α) − 1) = α,
hx, xi

(15)

where the inequality follows from Theorem 2.6 and noting that r still satisfies the theorem requirement. The result
follows by√adjusting the value of α after plugging this and equation (14) in equation (13), using the fact that kXk2 ≤
kXkF ≤ nkXk2 for any n × n matrix X.

5.2

Linear Regression

Our mechanism for L IN -R EG is presented below and the result is stated in Theorem 5.3. The main idea is to lift the
singular values of the input matrix above the threshold of Theorem 3.1(i).
p
I NITIALIZATION . On input parameters α, β, ε, δ, set r = O(d log(1/β)/α), s = 16r ln(2/δ)ε−1 ln(16r/δ), and
YA0 to be all zero matrix.
DATA - STRUCTURE UPDATE . On input a column c of an n × d matrix A at time epoch t, set the column vector
b :c = sec 0n+d A:c . Call PSG1 with input (A
b :c , r, n + d). Update the sketch of A by replacing the
A
column c of YAt−1 by the returned sketch to get the sketch YAt .

20


b i = 0d 0n+d bi . Call PSG1
A NSWERING QUERIES . On being queried with a vector bi , set the column vector b
b i , r, n + d) to get the sketch Yb . Compute a vector xi satisfying minx kYA xi − Yb k.
with input (b
i
t
i
Theorem 5.3. Let Ω be r × 2(n + d) matrix used by PSG, where r = O(d log(1/β)/α). Then the data-structure
generated above requires O(d2 α−1 κ log(1/β))
√ bits and allows to solve (α, β, τ )-L IN -R EG problem in an (ε, δ)differentially private manner with τ ≤ O(s2 nα).
Proof. The formal proof of Theorem 5.3 is identical to [36], which was refined in [10], modulo the analysis to consider
the lift of singular value and using the bound of [27] (Lemma 5.5). We first start with the intuition. For this, let us
b and B.
b Casted with these matrices, the result of [10] is as follows:
recall the main result of [10] with matrices A
b and B
b be the matrices as constructed in the mechanism stated in SecTheorem 5.4. Let Ω be as in Section 5.2, A
bX
b − Bk
b ≤ (1 + α)kA
bX
e − Bk,
b where
tion 5.2. Then with probability at least 1 − β, kA
b = arg min kΩ(AX
b − B)k
b 2
X
F

e = arg min k(AX
b − B)k
b 2.
X
F

X

T

X

T

b and B
b have the same block matrices
b T = I B . In other words, both A
b T = I A and B
Now note that A
e
b
b
for the first few rows. Since, X minimizes the value of (AX − B), it has entries 1 except for the last n positions.
b and B
b and the
Therefore, one way to look at the bound is that the multiplicative error is due to the last n entries of A
additive error is due to the approximation given by the sketch on the rest of the entries.
We follow the approach of [10], which is the refinement of [36]. Their approach to prove the utility bound for
b X
e − X)k
b F , where X
b is the right
linear regression works in two stages and we reiterate it here. We first bound kUT A(
e is the value of regression for the sketch of A
b and B,
b and U is left singular matrix of
value of the linear regression, X
b Let X
e 0 represents the last n entries of X
e (similarly define for X)
b and X
e 00 be the rest of the entries (similarly for
A.
b
X).
The following lemma given by [27] is key to our analysis. This is an equivalent for equation (15) with this value
of r. Note that a random Gaussian matrix of the order we considered is a Johnson-Lindenstrauss transform.
Lemma 5.5. Given r = O(d log(1/β)/α). Let U be any unitary matrix. If Ω satisfies the Johnson-Lindenstrauss
bound, then with probability at least 1 − β, we have kUΩT ΩUT − Ik2 ≤ α.
b X̃ − X)k
b F follows from the triangle inequality. The remainder of the
Using Lemma 5.5, the bound on kUT A(
proof follows from the Pythagorus theorem and an observation that A and U have the same column-space. A structural
b to A (respectively, from A
b to B) is that
property from our construction that we repeatedly use for moving from A
b
b
except for the block matrix formed by the last n rows, the rest of the entries of the matrices A and B are identical.
The details follows.
b Since
Let B be the matrix formed by the set of queries {b1 , · · · , bm } and U be the left singular matrix of A.
T
T
the columns of U is a set of orthonormal vectors, we have UU U = U and kU UCkF = kUCkF for any matrix
b X
b − X)k.
e
C. Therefore, it suffices for the utility bound to prove a bound on kUT A(
For this, we first prove that
b X
b − X)
e has a small norm. We have
UT ΩT ΩA(
b X
b − X)
e = UT ΩT ΩA(
b X
b − X)
e + UT ΩT Ω(B
b −A
b X)
b
UT ΩT ΩA(
b −A
b X).
e
= UT ΩT Ω(B
p
bX
b − B)
b = A
b T ΩT Ω(A
bX
b − B)
b = 0. Therefore, from Theorem 5.1 with α0 = α/d
This is because UT ΩT Ω(A
(since we chose r in Theorem 5.1 which differs by a factor of α and 1/d with respect to that in Theorem 5.3), we have
b X
b − X)k
e F = kUT ΩT Ω(B
b −A
b X)k
e F
kUT ΩT ΩA(
√
√
e 0 kF + τ
≤ αkB − AX
From the sub-additivity of the norm and property of conforming matrices (Lemma 2.1), we have
b X
b − X)k
e F ≤ kUT ΩT ΩA(
b X
b − X)k
e F + kUT ΩT ΩA(
b X
b − X)
e − UT A(
b X
b − X)k
e F
kUT A(
√
√
e 0 kF + τ + kUT ΩT ΩU − Ik2 · kUT A(
b X
b − X)k
e F
≤ αkB − AX

21

b X
b − X)k
e F ≤ 2√αkB − AX
b 0 kF + √τ . The utility
Using Lemma 5.5 and rearranging the terms, we get kUT A(
b and U are the same, X
e 0 is the optimal solution of the
proof is now immediate by observing that the column-space of A
T b e
T b e
e
b
b
b = 0), and the Pythagorus
regression (which from the normal equations of X implies U (AX − B) = A (AX − B)
theorem on the norms. More concretely, with probability at least 1 − 2β
b 0 − Bk2 = kAX
e 0 − Bk2 + kA(X
b0 − X
e 0 )k2
kAX
F
F
F
e 0 − BkF + τ.
≤ (1 + 4α)kAX
Adjusting and renaming the values of α and β, we get the claim of the theorem.

L IN -R EG IN OTHER MODELS . L IN -R EG has been also studied in the local privacy model by Duchi et al. [11] and in
the online learning model by Jain et al. [25] and Thakurata and Smith [38]. These models are different from ours and,
therefore, our results are incomparable to theirs. However, if one wishes to make any comparison with our result, we
give a brief overview of these results. The readers are welcome to compare
it with √
Theorem
 5.3. Jain et al. [25] used
√
R queries and T training data set to give a bound of Õ (R6 log(1/δ) n log1.5 T )/ εα3 . Thakurata and Smith [38]
gave a generic algorithm of which Jain et al. [25] is a special case. Duchi et al. [11] gave a characterization of L IN -R EG
in the local-privacy model.

5.3

Applications in Learning Theory

Due to lack of space, we do not cover few applications and just make a short note here. For example, PSG can be
used in manifold learning. Here, we can consider the streaming version of Baranuik and Wakin [2, Sec 3], where
the sampled points are streamed. We first perform a transformation as done in L IN -R EG on their sampled points to
boost the singular values of the matrix formed by the set of all sampled points and then use Gaussian matrix as the
Johnson-Lindenstrauss transform in their algorithm. The more formal treatment of this is done by Upadhyay [40].
Also, our mechanism for LRA can be easily compiled to give differentially private principal component analysis using
standard algorithms that use LRA on the input matrix as the first step. This achieves almost the same memory bound
as achieved by [30].

22

References
[1] Noga Alon, Yossi Matias, and Mario Szegedy. The space complexity of approximating the frequency moments.
In Proceedings of the twenty-eighth annual ACM symposium on Theory of computing, pages 20–29. ACM, 1996.
[2] Richard G Baraniuk and Michael B Wakin. Random projections of smooth manifolds. Foundations of computational mathematics, 9(1):51–77, 2009.
[3] Amos Beimel, Kobbi Nissim, and Eran Omri. Distributed private data analysis: Simultaneously solving how and
what. In Advances in Cryptology–CRYPTO 2008, pages 451–468. Springer, 2008.
[4] Jeremiah Blocki, Avrim Blum, Anupam Datta, and Or Sheffet. The Johnson-Lindenstrauss Transform Itself
Preserves Differential Privacy. In FOCS, pages 410–419, 2012.
[5] Jeremiah Blocki, Avrim Blum, Anupam Datta, and Or Sheffet. Differentially private data analysis of social
networks via restricted sensitivity. In ITCS, pages 87–96, 2013.
[6] Avrim Blum, Cynthia Dwork, Frank McSherry, and Kobbi Nissim. Practical privacy: the sulq framework. In
Chen Li, editor, PODS, pages 128–138. ACM, 2005.
[7] E Bura and R Pfeiffer. On the distribution of the left singular vectors of a random matrix and its applications.
Statistics &amp; Probability Letters, 78(15):2275–2280, 2008.
[8] Lynn E Cannon. A cellular computer to implement the kalman filter algorithm. Technical report, DTIC Document, 1969.
[9] Kamalika Chaudhuri, Anand D Sarwate, and Kaushik Sinha. Near-optimal differentially private principal components. In NIPS, pages 998–1006, 2012.
[10] Kenneth L. Clarkson and David P. Woodruff. Numerical linear algebra in the streaming model. In STOC, pages
205–214, 2009.
[11] John C Duchi, Michael I Jordan, and Martin J Wainwright. Local privacy and statistical minimax rates. In
Foundations of Computer Science (FOCS), pages 429–438. IEEE, 2013.
[12] Cynthia Dwork, Frank McSherry, Kobbi Nissim, and Adam Smith. Calibrating Noise to Sensitivity in Private
Data Analysis. In Shai Halevi and Tal Rabin, editors, TCC, volume 3876 of Lecture Notes in Computer Science,
pages 265–284. Springer, 2006.
[13] Cynthia Dwork, Moni Naor, Toniann Pitassi, and Guy N. Rothblum. Differential privacy under continual observation. In STOC, pages 715–724, 2010.
[14] Cynthia Dwork, Moni Naor, Toniann Pitassi, Guy N Rothblum, and Sergey Yekhanin. Pan-Private Streaming
Algorithms. In ICS, pages 66–80, 2010.
[15] Cynthia Dwork and Aaron Roth. The algorithmic foundations of differential privacy. Theoretical Computer
Science, 9(3-4):211–407, 2013.
[16] Cynthia Dwork, Guy N. Rothblum, and Salil P. Vadhan. Boosting and Differential Privacy. In FOCS, pages
51–60. IEEE Computer Society, 2010.
[17] Cynthia Dwork, Kunal Talwar, Abhradeep Thakurta, and Li Zhang. Analyze Gauss: Optimal Bounds for PrivacyPreserving Principal Component Analysis. In STOC, pages 11–20, 2014.
[18] Carl Eckart and Gale Young. The approximation of one matrix by another of lower rank. Psychometrika,
1(3):211–218, 1936.
[19] Arjun K Gupta and Daya K Nagar. Matrix variate distributions, volume 104. CRC Press, 1999.
[20] Nathan Halko, Per-Gunnar Martinsson, and Joel A Tropp. Finding structure with randomness: Probabilistic
algorithms for constructing approximate matrix decompositions. SIAM review, 53(2):217–288, 2011.
[21] Moritz Hardt and Eric Price. The noisy power method: A meta algorithm with applications. In Advances in
Neural Information Processing Systems, pages 2861–2869, 2014.

23

[22] Moritz Hardt and Aaron Roth. Beating randomized response on incoherent matrices. In STOC, pages 1255–1268,
2012.
[23] Moritz Hardt and Aaron Roth. Beyond worst-case analysis in private singular vector computation. In STOC,
pages 331–340, 2013.
[24] Piotr Indyk and Rajeev Motwani. Approximate nearest neighbors: Towards removing the curse of dimensionality.
In Proceedings of the Thirtieth Annual ACM Symposium on the Theory of Computing, Dallas, Texas, USA, May
23-26, 1998, pages 604–613, 1998.
[25] Prateek Jain, Pravesh Kothari, and Abhradeep Thakurta. Differentially Private Online Learning. In Shie Mannor,
Nathan Srebro, and Robert C. Williamson, editors, COLT, volume 23 of JMLR Proceedings, pages 24.1–24.34.
JMLR.org, 2012.
[26] Adam Kalai and Santosh Vempala. Efficient algorithms for online decision problems. Journal of Computer and
System Sciences, 71(3):291–307, 2005.
[27] Daniel M. Kane and Jelani Nelson. Sparser Johnson-Lindenstrauss Transforms. J. ACM, 61(1):4, 2014.
[28] Michael Kapralov and Kunal Talwar. On differentially private low rank approximation. In SODA, volume 5,
page 1. SIAM, 2013.
[29] Leon Mirsky. Symmetric gauge functions and unitarily invariant norms. The quarterly journal of mathematics,
11(1):50–59, 1960.
[30] Ioannis Mitliagkas, Constantine Caramanis, and Prateek Jain. Memory limited, streaming pca. In Advances in
Neural Information Processing Systems, pages 2886–2894, 2013.
[31] Elchanan Mossel, Ryan O’Donnell, and Rocco P Servedio. Learning juntas. In Proceedings of the thirty-fifth
annual ACM symposium on Theory of computing, pages 206–212. ACM, 2003.
[32] Robb J Muirhead. Aspects of multivariate statistical theory, volume 197. John Wiley &amp; Sons, 2009.
[33] C Radhakrishna Rao. Linear statistical inference and its applications, volume 22. John Wiley &amp; Sons,
2009.
[34] Mark Rudelson and Roman Vershynin. Smallest singular value of a random rectangular matrix. Communications
on Pure and Applied Mathematics, 62(12):1707–1739, 2009.
[35] Youcef Saad and Martin H Schultz. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM Journal on scientific and statistical computing, 7(3):856–869, 1986.
[36] Tamas Sarlos. Improved approximation algorithms for large matrices via random projections. In Foundations of
Computer Science, 2006. FOCS’06. 47th Annual IEEE Symposium on, pages 143–152. IEEE, 2006.
[37] Terence Tao. Topics in random matrix theory, volume 132. American Mathematical Soc., 2012.
[38] Abhradeep Guha Thakurta and Adam Smith. (nearly) optimal algorithms for private online learning in fullinformation and bandit settings. In Advances in Neural Information Processing Systems, pages 2733–2741,
2013.
[39] Jalaj Upadhyay. Random Projections, Graph Sparsification, and Differential Privacy. In ASIACRYPT (1), pages
276–295, 2013.
[40] Jalaj Upadhyay. Randomness Efficient Fast-Johnson-Lindenstrauss Transform with Applications in Differential
Privacy and Compressed Sensing . CoRR, abs/1410.2470, 2014.
[41] Thomas J Walsh, István Szita, Carlos Diuk, and Michael L Littman. Exploring compact reinforcement-learning
representations with linear regression. In Proceedings of the Twenty-Fifth Conference on Uncertainty in Artificial
Intelligence, pages 591–598. AUAI Press, 2009.
[42] Christopher Williams and Matthias Seeger. The effect of the input density distribution on kernel-based classifiers.
In Proceedings of the 17th International Conference on Machine Learning, number EPFL-CONF-161323, pages
1159–1166, 2000.

24

