Principal Component Analysis in the Local Differential Privacy Model

1 Introduction

Principal Component Analysis (PCA) is a fundamental technique for dimension reduction in statistics, machine learning, and signal processing. As of today, it remains as one of the most commonly used tools in applications, especially in social sciences [Costello and Osborne, 2005], financial econometrics [Aït-Sahalia and Xiu, 2017], medicine [Barber et al., 1975], and genomics [Lu and Xu, 2013].

With the rapid development of information technologies, big data now ubiquitously exist in our daily life, which need to be analyzed (or learned) statistically by methods like regression and PCA. However, due to the presence of sensitive data (especially those in social science, biomedicine and genomics) and their distributed nature, such data are extremely difficult to aggregate and learn from. Consider a case where health records are scattered across multiple hospitals (or even countries), it is challenging to process the whole dataset in a central server due to privacy and ownership concerns. A better solution is to use some differentially private mechanisms to conduct the aggregation and learning tasks. Differential Privacy (DP) (Dwork et al., 2006b) is a commonly-accepted criterion that provides provable protection against identification and is resilient to arbitrary auxiliary information that might be available to attackers.

<The research of this work was supported in part by NSF through grant CCF-1716400.

£Contact Author

Currently, there are mainly two user models available for differential privacy: the central model and the local one. In the central model, data are managed by a trusted central entity which is responsible for collecting them and for deciding which differentially private data analysis to perform and to release. A classical application of this model is the one of collecting census data. In the local model instead, each individual manages his/her proper data and discloses them to a server through some differentially private mechanisms. The server collects the (now private) data of each individual and combines them into a resulting data analysis. A classical example of this model is the one aiming at collecting statistics from user devices like in the case of Google's Chrome browser [Erlingsson et al., 2014], and Apple's iOS-10 [Near, 2018].

In the local model, two basic types of protocols are often used: interactive and non-interactive. [Smith et al., 2017] have recently investigated the power of non-interactive differentially private protocols. This type of protocols is more natural for the classical use cases of the local model: both projects from Google and Apple use the non-interactive model. Moreover, implementing efficient interactive protocols in such applications is more difficult due to the latency of the network. Despite being used in industry, the local model has been much less studied than the central one. Part of the reason for this is that there are intrinsic limitations in what one can do in the local model. As a consequence, many basic questions, that are well studied in the central model, have not been completely understood in the local model, yet.

In this paper, we study PCA under the non-interactive local differential privacy model and aim to answer the following main question.

What are the limitations and the (near) optimal algorithms of PCA under the non-interactive local differential privacy model?

We summarize our main contributions as follows:

1. We first study the k-subspace PCA problem in the low dimensional setting and show that the minimax risk (measured by the squared subspace distance) under � non-interactive local differential privacy (LDP) is lower bounded by�. �1�k+1pk .�k*�k+1/2n�2 /, where p is the dimensionality of the data and n is the number of data records. Moreover, we prove that the term�. pk n�2 / is optimal.
2. An undesirable issue of the above result is that the er-

ror bound could be too large in high dimensions (i.e., p , n). In such scenarios, a natural approach is to impose some additional structural constraints on the leading eigenvectors. A commonly used constraint is to assume that the leading eigenvectors are row sparse, which is refereed as sparse PCA in the literature and has been studied intensively in recent years [Vu et al., 2013b; Cai et al., 2013; Vu et al., 2013a]. Thus, for the high dimensional case, we consider the sparse PCA under the noninteractive local model and show that the private minimax risk (measured by the squared subspace distance) is lower bounded by �. �1�k+1 .�k*�k+1/2 ks log p n�2 /, where �1; �k and �k+1 are eigenvalues and s is the sparsity parameter of the eigenvectors. We also give an algorithm to achieve a near optimal upper bound of O� �21 .�k*�k+1/2 s2 logp n�2 �.

2 RelatedWork

There is a vast number of papers studying PCA under differential privacy, starting from the SULQ framework [Blum et al., 2005], [Dwork et al., 2014; Chaudhuri et al., 2013; Jiang et al., 2016; Gonem and Gilad-Bachrach, 2018; Ge et al., 2018; Balcan et al., 2016]. We compare only those private PCA results in distributed settings.

For the low dimensional case, Balcan et al. [Balcan et al., 2016] studied the private PCA problem under the interactive local differential privacy model and introduced an approach based on the noisy power method. They showed an upper bound which is suitable for general settings, while ours is mainly for statistical settings. It is worth pointing out that the output in [Balcan et al., 2016] is only an O.k/-dimensional subspace, instead of an exact k-dimensional subspace; thus their result is incomparable with ours. Moreover, we provide, in this paper, a lower bound on the � non-interactive private minimax risk.

For the private high dimensional sparse PCA, the work most closely related to ours is the one by Ge et al. [Ge et al., 2018]. The authors in this paper proposed a noisy iterative hard thresholding power method, which is an interactive LDP algorithm and proved an upper bound of O� �1�k .�k*�k+1/2 s.k+log p/ n.1*� 14 / � for their method, where � is a param-

eter related to �. Specifically, they showed that there exists some 'Privacy Free Region'. However, several things need to be pointed out. Firstly, our method is for general � ¸ .0; 1] and non-interactive settings, while Ge et al. considered the interactive setting with more restricted �. Secondly, the assumptions in our paper are less strict than the ones in [Ge et al., 2018]. Finally, we provide a lower bound on the private minimax risk.

The optimal procedure in our paper is based on perturbing the covariance by Gaussian matrices, which has been studied in [Dwork et al., 2014]. However, there are some major differences; firstly, we show the optimality of our algorithm under the non-interactive local model using subspace distance as the measurement, while [Dwork et al., 2014] showed the optimality under the .�; �/ central model using variance as the measurement. It is notable that in [Dwork et al., 2014] the authors also provided an upper bound on the subspace distance. However, the lower bound is still unknown. Secondly, while the optimal algorithm for the low dimensional case is quite similar, we extend it to the high dimensional case. The optimal procedure in the high dimensional sparse case is quite different from that in [Dwork et al., 2014]. Thirdly, in this paper,since we focus the statistical setting while [Dwork et al., 2014] considered the general setting, the upper bound results are incomparable.

3 Preliminaries

3.1 Classical Minimax Risk

Since all of our lower bounds are in the form of private minimax risk, we first introduce the classical statistical minimax risk before discussing the locally private version.

Let P be a class of distributions over a data universe X. For each distribution P ¸ P, there is a deterministic function �.P / ¸ �, where � is the parameter space. Let � : ��� :› R+ be a semi-metric function on the space � and � : R+ › R+ be a non-decreasing function with �.0/ = 0 (in this paper, we assume that � is the subspace distance and �.x/ = x2 unless specified otherwise).

We further assume that ^Xi'ni=1 are n i.i.d observations drawn according to some distribution P ¸ P, and �� : X n › � is some estimator. Then, the minimax risk in metric ��� is defined by the following saddle point problem:

where the supremum is taken over distributions P ¸ P and the infimum over all estimators ��.

3.2 Local Differential Privacy and Private Minimax Risk

Since we consider the non-interactive local model in this paper, we will follow the definitions in [Duchi et al., 2013].

We assume that ^Zi'ni=1 are the private observations transformed from ^Xi'ni=1 through some privacy mechanisms. When Zi depends only on Xi, the mechanism is called noninteractive and in this case we have a simpler form for the conditional distributions Qi.Zi � Xi = xi/. We now define local differential privacy by restricting the conditional distribution.

Definition 1 ([Duchi et al., 2013]). For a given privacy parameter �; � > 0, we say that the random variable Zi is an .�; �/ non-interactively locally differentially private (LDP) view of Xi if for any xi; x¤i ¸ X:

When � = 0, we call it � non-interactively LDP view. We say that the privacy mechanismQ = ^Qi'ni=1 is .�; �/ (�/ noninteractively locally differentially private (LDP) if each Zi is an .�; �/ (�/ non-interactively LDP view.

For a given privacy parameter � > 0, let Q� be the set of conditional distributions that have the �-LDP property. For a given set of samples ^Xi'ni=1 , let ^Zi'ni=1 be the set of observations produced by any distribution Q ¸ Q�. Then, our estimator will be based on ^Zi'ni=1 , that is, ��.Z1;5;Zn/. This yields a modified version of the minimax risk:

From the above definition, it is natural to seek the mechanism Q ¸ Q� that has the smallest value for the minimax risk. This allows us to define functions that characterize the optimal rate of estimation in terms of privacy parameter �.

Definition 2. Given a family of distributions �.P/ and a privacy parameter � > 0, the � non-interactive private minimax risk in the metric ��� is:

where Q� is the set of all � non-interactively locally differentially private mechanisms.

3.3 Locally Private k-dimensional PCA

Let X ¸ Rp a random vector with mean 0 and covariance matrix �. k-dimensional PCA is to find a k dimensional subspace that optimizes the following problem:

where Gp;k is the Grassmann manifold of k-dimensional subspaces of Rp , and �G is the projection of G. There always exists at least one solution; consider � = ‡pj=1 �jvjvTj , where �1 g �2 g 5; �p g 0 are the eigenvalues of � and v1; v2;5; vp ¸ Rp are the corresponding eigenvectors. If �k g �k+1, then the k-dimensional principal subspace of �, i.e. the subspace S spanned by v1;5; vk solves the above optimization problem, where the orthogonal projector of S is given by �S = VkV T k , where Vk = [v1;5; vk] ¸ Vp;k, Vp;k is the set of all p � k orthogonal matrices. For simplicity we denote S = col.Vk/, where col.M/ denotes the subspace spanned by the columns vectors of M.

In practice, � is unknown, and the only thing that we have is the set of observation data records ^X1;5;Xn', which are i.i.d sampled from X. Thus, the problem of (noninteractively) locally differentially private PCA is to find a kdimensional subspace S priv which is close to S, where the algorithm that outputs S priv must be � (non-interactively) locally differentially private. After obtaining a private estimator S priv , there are multiple ways to measure the success, such as variance guarantee [Dwork et al., 2014], low rank approximation error [Kapralov and Talwar, 2013], etc. In this paper, we will use the subspace distance as the measurement [Dwork et al., 2014; Ge et al., 2018].

Let S and S ¤ be two k-dimensional subspaces in Rp . Also denote by E and F, respectively, the orthogonal matrix corresponds to S and S ¤ . That is, E = V V T and F = WW T for some orthogonal matrices V ¸ Vp;k and W ¸ Vp;k. Then, the squared subspace distance between S and S ¤ is defined by the following [Stewart, 1990]:

where æ � æF is the Frobenious norm. For simplicity, we will overload notation and write sin �.S; S ¤ / = sin �.V ; W /.

4 LowDimensional Case

In this section, we focus on the general case and always assume n g p. We first derive a lower bound of the � noninteractive private minimax risk using the squared subspace distance as the measurement. By the definition of the �- private minimax risk, it is important to select an appropriate class of distributions.

4.1 Class of Distributions

We assume that the observations ^Xi'ni=1 satisfies: X = � 1 2 Z; where Z ¸ Rp is a random vector satisfying equations EZ = 0 and Var.Z/ = Ip. Also, Z is sub-Gaussian and æZæ 2 f 1 where

which means that all the one-dimensional marginals of Xi have sub-Gaussian tails. We need to note that this assumption on X is commonly used in many papers on PCA in statistical settings, such as [Vu et al., 2013b; Ge et al., 2018].

In the study of private PCA, it is always assumed that the l2 norm of each Xi is bounded by 1, as in [Dwork et al., 2014][Ge et al., 2018]. For convenience, we relax this assumption in the following way; for the random vector X ¸ Rp , we assume that æXæ2 f 1 with a probability at least 1*e*�.p/ .

Next, we give assumptions on the population covariance matrix �. Firstly, we assume that for the target k-dimensional subspace, �k *�k+1 > 0 so that the principal subspace is well defined. Next, we define the effective noise variance, which is proposed in [Vu et al., 2013b] and [Cai et al., 2013]:

For a given constant �2 > 0, we assume that �2 k f �2 . Also, we denote the collection of distributions by P.k; �2 /.

4.2 Main Results

The next theorem shows a lower bound of � non-interactive private minimax risk under squared subspace distance.

Theorem 1. Let ^Xi'ni=1 be samples from P ¸ P.k; �2 /. If p4 f k f 3p 4 , � ¸ .0; 12 ] and n g �� 1 �2 �1�k+1 .�k*�k+1/2 min^k; p * k'�, then the � non-interactive private minimax risk in the metric of squared subspace distance satisfies:

Remark 1. We note that for the non-private case, the minimax risk is lower bounded by �� �1�k+1 .�k*�k+1/2 kp n � [Cai et al., 2013]. Thus, in this case, the impact of the local differential privacy is to change the number of efficient sample from n to n�2 . However, the collection of the considered distributions needs another assumption which says that æXæ2 is bounded by 1 with high probability. This is not necessary in the nonprivate case [Cai et al., 2013], but needed in ours for showing the upper bound.

We also note that in the central differential privacy model, [Dwork et al., 2014] showed that the lower bound of the kdimensional PCA is ��. kp log. 1 � / n2�2 / for .�; �/-differential privacy. However, this lower bound is measured by the variance of X = .XT 1 ;XT 2 ;5;XT n /T ¸ Rn�p , not the squared subspace distance used in this paper. Although [Dwork et al., 2014] gave an upper bound of O. kp log.1_�/ .�2k*�2k+1 /n2�2 / in the general setting using the squared subspace distance as measurement, it is still unknown whether the bound is optimal. Also, their lower bound omits the parameters related to the eigenvalues. For the � differential privacy in the central model, [Chaudhuri et al., 2013] showed that the lower bound is �. p2 n2�2.�1*�2/2 / in the special case of k = 1. However, it is still unknown for the general case of k. Thus, from the above discussion, we can see that the lower bound of � non-interactively locally differentially private PCA is similar to the .�; �/ differentially private PCA in the central model.

One of the main questions is whether the lower bound in Theorem 1 is tight. In the following, we show that the term �. pk n�2 / is tight. By our definition of the parameter space, we know that for any X � P ¸ P.�2 ; k/, æXæ2 f 1 with high probability. Thus, we always assume that the event of each æXiæ2 f 1 holds. Note that this assumption also appears in [Ge et al., 2018; Dwork et al., 2014; Balcan et al., 2016]. The idea is the same as in [Dwork et al., 2014], where each Xi perturbs its covariance and aggregates the noisy version of covariance, see Algorithm 1 for details.

Theorem 2. For any �; � > 0, Algorithm 1 is .�; �/ (noninteractively) locally differentially private. Furthermore, with probability at least 1 * e*C1 * 1 pC2 , the output satisfies:

where C1; C2 are some universal constants.

Algorithm 1 Local Gaussian Mechanism

Input: data records ^Xi'ni=1 � P n for P ¸ P.�2 ; k/, and for i ¸ [n]; æXæ2 f 1. �; � are the privacy parameters.

• 1: for Each i ¸ [n] do
• 2: Denote �Xi �XT i = XiXT i + Zi, where Zi ¸ Rp�p is a symmetric matrix where the upper triangle ( including the diagonal) is i.i.d samples fromN.0; �2 1 /; here �2 1 = 2 ln.1:25_�/ �2 , and each lower triangle entry is copied from its upper triangle counterpart.
• 3: end for
• 4: Compute �S = 1n ‡ni=1 �Xi �XT i .
• 5: Output col. � Vk/ where � Vk ¸ Rp�k is the principal rank k subspace of �S.

In Theorem 7 of [Dwork et al., 2014], the authors provided a similar upper bound for the .�; �/-differential privacy in the central model. However, they need to assume that the eigenvalues satisfy the condition �2k * �2k+1 = !.øp/, which is not needed in our Theorem 2.

From Theorems 1 and 2, we can see that there is still a gap of O. �1 �k+1 / between the lower and upper bounds. We leave it as an open problem to determine whether these bounds are tight or not.

5 High Dimensional Sparse Case

From Theorem 1, we can see that for the high dimensional case, i.e. p , n, the bound in (2) becomes trivial. Thus, to avoid this issue, we need some additional assumption on the parameter space. One of the commonly used assumption is sparsity. There are many definitions of sparsity on PCA and we use the row sparsity in this paper, which has also been studied in [Vu et al., 2013b; Cai et al., 2013; Ge et al., 2018].

We first define the .2; 0/-norm of a p � k matrix A as the usual l2 norm of the vector of row-wise l0 norms of A:

where aj< denotes the j-th row of A. Note that æ � æ2;0 is coordinate independent, i.e. æAOæ2;0 = æAæ2;0 for any orthogonal matrix O ¸ Rk�k . We define the row sparse space as follows.

Definition 3. Let s be the sparsity level parameter satisfying the condition of k f s f p. The s-(row) sparse subspace is defined as follows

We define our parameter space, P.s; k; �2 /, to be the same as in the previous section with an additional condition that S ¸ M0.s/, where S is the k-dimensional principal subspace of covariance matrix �.

Below, we will first derive a lower bound of the noninteractive locally differentially private PCA in the high dimensional sparse case.

Theorem 3. Let ^Xi'ni=1 be the observations sampled from a distribution P ¸ P.s; k; �2 /. If the privacy parameter � ¸ .0; 1 2 ], n g C.s * k/ �2 .k+log p/ �2 for a universal constant C > 0, and 2k f s * k f p * k, then for all k ¸ [p] satisfying the condition of k f p * 4, the � non-interactive private minimax risk in the metric of squared subspace distance satisfies the following

Note that in the non-private case, the optimal minimax risk is O� �1�k+1 .�k*�k+1/2 s.k+log p/ n �. Thus, same as in the low dimensional case, the impact of the privacy constraint is to change the efficient samples from n to n�2 .

Next, we consider the upper bound. In the non-private case, the optimal procedure is to solve the following NP-hard optimization problem [Vu et al., 2013b]:

where S is the empirical covariance matrix. Our upper bound is based on (4). However, instead of solving (4) on the perturbed version of the empirical covariance matrix, we perturb the covariance matrix and solve the following optimization problem on the convex hull of the constraints in (4), that is:

where < S; X >= Tr.SXT /. Note that the constraints in (5), which is called Fantope [Bhatia, 2013][Vu et al., 2013a], is the convex hull of the constrains in (4). Also, since the constraints in (5) only guarantees that the rank of the output is g k, the output �X needs not to be a matrix with exact rank of k. Thus, in order to obtain a proper k-dimensional subspace, we just output the k-PCA of �X.

Algorithm 2 Local Gaussian Mechanism-High Dimension

Input: data records ^Xi'ni=1 � P n for P ¸ P.s; �2 ; k/, and for i ¸ [n]; æXæ2 f 1. �; � are privacy parameters. � > is a constant.

• 1: for Each i ¸ [n] do
• 2: Denote �Xi �XT i = XiXT i + Zi, where Zi ¸ Rp�p is a symmetric matrix where the upper triangle ( including the diagonal) is i.i.d samples fromN.0; �2 1 /; here �2 1 = 2 ln.1:25_�/ �2 , and each lower triangle entry is copied from its upper triangle counterpart.
• 3: end for
• 4: Compute �S = 1n ‡ni=1 �Xi �XT i .
• 5: Get the optimal solution �X in (5) or do as the followings 6: Setting Y .0/ = 0;U .0/ = 0
• 7: for t = 1; 2;5do

• 9: Y .t+1/ = S�_�.X.t+1/ + U .t/ / where S is the entrywise soft thresholding operator defined as S�_�.x/ = sign.x/ max.�x� * �_�; 0/.
• 10: U .t+1/ = U .t/ +X.t+1/ * Y .t+1/
• 11: Return Y .t/
• 12: end for
• 13: Let k-dimensional principal component of �X or Y .t/ be � Vk, output �S = col. � Vk/.

Theorem 4. For any given 0 < �; � < 1, if ^Xi'ni=1 � P n for P ¸ P.s; �2 ; k/ and æXæ2 f 1 for all i ¸ [n], then the solution to the optimization problem (5) is .�; �/ non-interactive locally differentially private. Moreover, if let � Vk denote the k-dimensional principal component subspace of �X and set � f O.�1t log p n�2 /, then with probability at least 1 * 2 p2 * 1 pc , the following holds

where c is a universal constant.

Since the optimization problem (5) is convex, we can follow the approach in [Vu et al., 2013a] to solve it by using ADMM method (see Algorithm 2 for the details).

Comparing with the lower bound of the private minimax risk in Theorem 3, we can see that the bound in Theorem 4 is roughly larger than the optimal rate by a factor of O. �1 �k+1 sk /.

This means that the upper bound is only near optimal [Vu et al., 2013a]. A remaining open problem is to determine whether it is possible to get a tighter upper bound that does not contain the term of sk in the gap.

6 Experiments

6.1 LowDimensional Case

For synthetic datasets, we generate the data samples ^Xi'ni=1 independently from a multivariate Gaussian distribution N.0; �/, where � = � 5p.�+1// V V T + 1 5p.�+1/ Ip for V ¸ Vp;k. It can be shown that æXiæ2 f 1 ¯i ¸ [n] with high probability (see supplemental material). We choose n = 105 , p = 40, k = ^5; 10; 15; 20', � = 0:5, � = 10*4 , and � = 1. For real world datasets, we run Algorithm 1 on Covertype and Buzz datasets [Dheeru and Karra Taniskidou, 2017] with normalized rows for each dataset. The error is measured by the subspace distance æ � Vk � V T k * VkV T k æF . For each experiment, we repeat 20 times and take the average as the final result.

Figure 1 is the result for the synthetic datasets. Figure indicates that 1) the error deceases as the sample size increases or � increases (i.e., becomes less private); 2) the error increases as the dimensionality p increases or the dimensionality k of the target subspace increases.

6.2 High Dimensional Case

For the high dimensional case, we consider the same distributions as in the low dimensional case and generate the target subspace V in the following way. For a given sparsity parameter s, we first generate a random orthogonal matrix � V ¸ Rs�k , then pad it with rows of zeros, and finally randomly permute the matrix. We set k = 10, n = 2000, p = 400, s = ^15; 20; 40; 80' and � = 1.

Besides the synthetic datasets, we also test our algorithm on some real world datasets in [Dheeru and Karra Taniskidou, 2017] and [Li et al., 2017]. We can see that 1) as the term of kn increases (n decreases), the error increases accordingly; 2) the error slightly increases when the dimensionality p increases, which is due to the fact that the upper bound in Theorem 4 depends only logarithmically on p (i.e., logp); 3) the error decreases when � increases. Table 1 and 2 show the results of the error with different sparsity and privacy, respectively.