High-Dimensional Differentially-Private EM Algorithm: Methods and Near-Optimal Statistical Guarantees

The EM algorithm is a widely used algorithm to compute the maximum likelihood estimator when there are latent or unobserved variables in the model. We first introduce the standard EM algorithm. Let $\bm{Y}$ and $\bm{Z}$ be random variables taking values in the sample spaces $\mathcal{Y}$ and $\mathcal{Z}$. For each pair of data $(\bm{Y},\bm{Z})$, we assume that only $\bm{Y}$ is observed, while $\bm{Z}$ remains unobserved. Suppose that the pair $(\bm{Y},\bm{Z})$ has a joint density function $f_{\bm{\beta}^{*}}(\bm{y},\bm{z})$, where $\bm{\beta}^{*}$ is the true parameter that we would like to estimate. Let $h_{\bm{\beta}^{*}}(\bm{y})$ be the marginal density function of the observed variable $\bm{Y}$. Then, we can write $h_{\bm{\beta}^{*}}(\bm{y})$ by integrating out $\bm{z}$

The goal of the EM algorithm is to obtain an estimator of $\bm{\beta}^{*}$ through maximizing the likelihood function. Specifically, suppose we have $n$ $i.i.d.$ observations of $\bm{Y}$: $\bm{y}_{1},\bm{y}_{2},...\bm{y}_{n}\sim h_{{\bm{\beta}}^{*}}(\bm{y})$, we aim to use EM algorithm to maximize the log-likelihood $l_{n}(\bm{\beta})=\sum_{i=1}^{n}\log h_{\bm{\beta}}(\bm{y}_{i})$, and get an estimator of ${\bm{\beta}}^{*}$. In many latent variable models, it is generally difficult to evaluate the log-likelihood $l_{n}(\bm{\beta})$ directly, but relatively easy to compute the log-likelihood for the joint distribution $f_{\bm{\beta}}(\bm{y},\bm{z})$. Such situations are in need of EM algorithms. Specifically, for a given parameter ${\bm{\beta}}$, let $k_{\bm{\beta}}(\bm{z}|\bm{y})$ be the conditional distribution of $\bm{Z}$ given the observed variable $\bm{Y}$, that is, $k_{\bm{\beta}}(\bm{z}|\bm{y})=\frac{f_{\bm{\beta}}(\bm{y},\bm{z})}{h_{\bm{\beta}}(\bm{y})}.$

The EM algorithm uses an iterative approach motivated by Jensen’s inequality. Suppose that in the $t$-th iteration, we have obtained ${\bm{\beta}}^{t}$ and would like to update it to ${\bm{\beta}}^{t+1}$ with a larger log-likelihood. The log-likelihood evaluated at $\bm{\beta}^{t+1}$ can always be lower bounded, as shown in the following expression.

with equality holds when ${\bm{\beta}}^{t+1}={\bm{\beta}}^{t}$. The basic idea of EM algorithm is to successively maximize the lower bound with respect to ${\bm{\beta}}^{t+1}$ in (2). In the E-step, we evaluate the lower bound in (2) at the current parameter ${\bm{\beta}}^{t}$. Since the second term in (2) only depends on $\bm{\beta}^{t}$, which is fixed given the current ${\bm{\beta}}^{t}$, we only need to evaluate $Q_{n}$ in the E-step. Then, in the M-Step, we calculate a new parameter ${\bm{\beta}}^{t+1}$ which moves towards the direction that maximizes $Q_{n}$, so the lower bound in (2) becomes larger when we update ${\bm{\beta}}^{t}$ to ${\bm{\beta}}^{t+1}$. We use the new parameter ${\bm{\beta}}^{t+1}$ at the $(t+1)$-th iteration and continue the E-step and M-step iteratively until convergence.

In the $t$-th iteration, the M-step in the standard EM algorithm maximizes $Q_{n}(\cdot;\bm{\beta}^{t})$ (Dempster et al., 1977), that is, $\bm{\beta}^{t+1}=\text{argmax}_{\bm{\beta}}Q_{n}(\bm{\beta};\bm{\beta}^{t})$. However, sometimes it is computationally expensive to compute the maximizer directly. As an alternative, the gradient EM (Balakrishnan et al., 2017) was proposed by taking one-step update of the gradient ascent $\bm{\beta}^{t+1}=\bm{\beta}^{t}+\eta\cdot\nabla Q_{n}(\bm{\beta}^{t};\bm{\beta}^{t})$ in the M-step. When $Q_{n}(\bm{\beta};\bm{\beta}^{\prime})$ is strongly convex with respect to ${\bm{\beta}}$, this gradient EM approach is shown to reach the same statistical guarantee as that of the standard EM approach (Balakrishnan et al., 2017). In the high-dimensional setting, when the data dimension is much larger than the sample size and the true parameter is sparse, people have come up with different variants of EM algorithms. For example, in the M-step, the maximization approach is generalized to the regularized maximization (Yi and Caramanis, 2015; Cai et al., 2019a; Zhang et al., 2020) and the gradient approach is generalized to the truncated gradient (Wang et al., 2015).

In this section, we introduce the concepts and properties of differential privacy. These properties will play an important role in the design of the DP EM algorithm. First, the formal definition of differential privacy is given below.

According to the definition, the two parameters $\epsilon,\delta$ control the privacy level. With smaller $\epsilon$ and $\delta$, the privacy constraint becomes more stringent. Intuitively speaking, this definition suggests that for a DP algorithm $M$, an adversary cannot distinguish if the original dataset is $\bm{X}$ or $\bm{X}^{\prime}$ when $\bm{X}$ and $\bm{X}^{\prime}$ is adjacent, implying that the information of each individual is protected after releasing $M$.

We then list several useful facts of designing DP algorithms. To create a DP algorithm, the arguably most common strategy is to add random noise to the output. Intuitively, the scale of the noise can not be too large, otherwise the accuracy of the output will be sacrificed. This scale is characterized by the sensitivity of the algorithm.

For algorithms with finite $\ell_{1}$-sensitivity or $\ell_{2}$-sensitivity, we can achieve differential privacy by adding Laplace noises or Gaussian noises respectively.

These two mechanisms are computationally efficient, and are typically used to build more complicated algorithms. In the following, we introduce the post-processing and composition properties of differential privacy, which enable us to design complicated DP algorithms by combining simpler ones.

In the following section, we will see that the above two properties are particularly useful in the design of the DP EM algorithm.

Suppose we have $i.i.d$ data sampled from the latent variable model (1) and aim to use the EM algorithm to find the maximum likelihood estimator in a DP manner. Since the EM algorithm is an iterative approach where the $t$-th iteration takes as input the ${\bm{\beta}}^{t-1}$ from the M-step in the last iteration, it suffices to make each ${\bm{\beta}}^{t}$ differentially private to ensure the privacy guarantee of the final output.

Our algorithm relies on two key designs in the M-step. First, we use the gradient EM approach, and in the gradient update stage, we introduce a truncation step on the gradient to control the sensitivity of the gradient update, and thus we can appropriately determine the scale of noise we need to achieve differential privacy. Secondly, we propose to apply the noisy hard-thresholding algorithm (NoisyHT) (Dwork et al., 2018) to pursue sparsity and achieve privacy at the same time. The NoisyHT algorithm is defined as follows,

In the last step, $P_{S}(\bm{u})$ denotes the operator that makes $\bm{u}_{S^{c}}=0$ while preserving $\bm{u}_{S}$. A great advantage of this algorithm is that when the vector $\bm{v}=v(\bm{X})$ has bounded $\ell_{\infty}$ sensitivity $\lambda$, the algorithm is guaranteed to be DP, as we can see in the lemma below.

Another important property of the NoisyHT algorithm is its accuracy. Specifically, for the coordinates that are not chosen by the NoisyHT algorithm, the $\ell_{2}$ norm of these coordinates is upper bounded by that of the coordinates with the same size but chosen by the NoisyHT algorithm up to some error term. The formal statement is summarized in the following with the proof given in Appendix C.1.

After introducing the NoisyHT algorithm, we now proceed to the development of the DP EM algorithm. We update the estimator $\hat{\bm{\beta}}$ based on the NoisyHT algorithm in each M-step after truncation. Such a modified M-step guarantees that the output ${\bm{\beta}}^{t}$ is sparse and differentially private in each iteration. We also note here that we use the sample-splitting in the algorithm, which makes ${\bm{\beta}}^{t}$ independent with the batch of samples in $t$-th iteration and helps control the sensitivity of the gradient. The algorithm is summarized below:

In the above algorithm, the truncation operator $f_{T}(\nabla Q_{n}(\bm{\beta};\bm{\beta}))$ with truncation level $T$ is defined as $f_{T}(\nabla Q_{n}(\bm{\beta};\bm{\beta}))=\frac{1}{n}\sum_{i=1}^{n}h_{T}(\nabla q_{i}(\bm{\beta};\bm{\beta}))$, where $q_{i}(\bm{\beta};\bm{\beta})=\int_{\mathcal{Z}}k_{{\bm{\beta}}}(\bm{z}_{i}|\bm{y}_{i})\log f_{\bm{\beta}}(\bm{y}_{i},\bm{z}_{i})d\bm{z}_{i}$, and $h_{T}$ denotes some generic truncation function with $\ell_{\infty}$ norm upper bounded by $T$. Later in the application to different specific models, we will specify the form of $h_{T}$ for each model. The next lemma shows Algorithm 2 is $(\epsilon,\delta)$-DP.

The proof of Lemma 3.3 is given in the Appendix C.2. In theory, we shall take the truncation level $T$ and number of iterations $N_{0}$ to be the order of $\sqrt{\log n}$ and $\log n$, respectively, as we will show in the next section. For the sparsity level $\hat{s}$, we should choose it to have the same order of the true sparsity level $s^{*}$. While the true parameter $s^{*}$ is unknown in practice, $\hat{s}$ could be chosen through cross-validation. Given the differential privacy guarantees, we are going to analyze the utility of this algorithm in the next section.

In this section, we analyze the theoretical properties of the proposed high-dimensional DP EM algorithm. Before we lay out the main results, we first introduce four technical conditions. The first three conditions are standard in the literature of EM algorithms; for example, see (Balakrishnan et al., 2017; Wang et al., 2015; Zhu et al., 2017; Wang et al., 2020). The fourth condition is needed to control the sensitivity in the design of DP algorithms. We will verify these four conditions in the specific models in the next section.

This condition implies that if $\bm{\beta}$ is close to the true parameter $\bm{\beta}^{*}$, then the gradients $\nabla Q(\bm{\beta};\bm{\beta}^{*})$ and $\nabla Q(\bm{\beta};\bm{\beta})$ should also be close, which implies gradient stability.

The Concavity-Smoothness condition indicates that when the second argument of $Q(\cdot;\cdot)$ is $\bm{\beta}^{*}$, $Q(\cdot;\cdot)$ is a well-behaved convex function that it can be upper bounded and lower bounded by two quadratic functions. This condition ensures the geometric decay of the optimization error in the statistical analysis.

In this condition, the statistical error is quantified by the $\ell_{\infty}$ norm between the population quantity $\nabla Q(\bm{\beta};\bm{\beta})$ and its sample version $\nabla Q_{n}(\bm{\beta};\bm{\beta})$. Such a bound is different from the $\ell_{2}$ norm bound considered in the classic low-dimensional DP EM algorithms (Wang et al., 2020). In the high-dimensional setting, although for each index, the statistical error is small, the $\ell_{2}$ norm can still be quite large. This fine-grained $\ell_{\infty}$ bound enables us to iteratively quantify the statistical accuracy when using the NoisyHT in the M-step.

The Truncation-Error condition quantifies the error caused by the truncation step in Algorithm 2. Intuitively, when $T$ is large, the truncation error $\xi$ can be very small, while leading to larger sensitivity and larger injected noise to ensure privacy. We need to carefully choose $T$ to strike a balance between the statistical accuracy and privacy cost. Below we show the main result of this section, with the detailed proof given in Section B.1.

To interpret this result, let us discuss the four terms in (3). The first term, $\kappa^{t/2}\cdot||{\bm{\beta}}^{0}-\bm{\beta}^{*}||_{2}=O(\kappa^{t/2}\cdot\|{\bm{\beta}}^{*}\|_{2})$, is the optimization error. With $\kappa\in(0,1)$, this term shrinks to zero at a geometric rate when the iteration number $t$ is sufficiently large. The second term is of order $\sqrt{s^{*}}\cdot\alpha$ when $\hat{s}$ is chosen as the same order of $s^{*}$, this is the statistical error caused by finite samples. We will further show that in some specific models, $\alpha$ is of the order $\sqrt{\log d\cdot\log n/n}$ and makes the second term to be $O(\sqrt{s^{*}\log d\cdot\log n/n})$. The third term $K\cdot\frac{s^{*}\log d\cdot\sqrt{\log(1/\delta)}\log^{3/2}n}{n\epsilon}$ can be seen as the cost of privacy, as this error is introduced by the additional requirement that the output needs to be $(\epsilon,\delta)$-DP. This term becomes larger when the privacy constraint becomes more stringent ($\epsilon,\delta$ become smaller). Typically, in practice, we choose $\delta=O(1/n^{a})$ for some $a\geq 1$ and $\epsilon$ a small constant. This implies that when $d,s^{*}$ and $n$ satisfy $\frac{s^{*}\log d\cdot(\log n)^{2}}{n}\cdot\frac{\log{1/\delta}}{\epsilon^{2}}=o(1)$, the cost of privacy will be negligible comparing to the statistical error. In this case, we can gain privacy for free in terms of convergence rate. The fourth term is due to the truncation of the gradient. Under the Truncation-Error condition with an appropriately chosen truncation parameter, this term reaches a convergence rate dominated by the statistical error up to logarithm factors.

In this subsection, we first apply the results in Section 3 to the Gaussian mixture model. By verifying the conditions in Theorem 3.8, we establish the convergence rate of the DP estimation in the Gaussian mixture model.

In a standard Gaussian mixture model, we assume:

where $\bm{Y}$ is a $d$-dimensional output and $\bm{e}\sim N(\bm{0},\sigma^{2}\bm{I}_{d})$. In this model, $\bm{\beta}^{*}$ and $-{\bm{\beta}}^{*}$ are the $d$-dimensional vectors representing the population means of each class, and $Z$ is a class indicator variable with ${\mathbb{P}}(Z=1)=1/2$ and ${\mathbb{P}}(Z=-1)=1/2$. Note that $Z$ is a hidden variable and independent of $\bm{e}$. In the high dimensional setting, we assume $\bm{\beta}^{*}$ to be sparse.

Let $\bm{y}_{1},\bm{y}_{2}...\bm{y}_{n}$ be $n$ $i.i.d$ samples from the Gaussian mixture model. Using the framework of the EM method introduced in Section 3, we need to calculate

where

Then, for the $M$-step in the $t$-th iteration given $\bm{\beta}^{t}$, the update rule is given by

Given this expression, we now present the DP estimation in the high-dimensional Gaussian mixture model by applying Algorithm 2. The algorithm is presented in Algorithm 3.

To derive the converge rate of $\hat{\bm{\beta}}$, we need to first verify Conditions 3.4-3.7. Conditions 3.4-3.6 are standard in the literature of EM algorithms (Balakrishnan et al., 2017; Wang et al., 2015). We adapt them into our setting and state the results altogether below.

The detailed proof of Lemma 4.1 is given in Appendix C.3. Given these verified conditions, the following theorem establishes the results for the DP estimation in the high-dimensional Gaussian mixture model.

The proof of Theorem 4.2 is given in Section B.2. Similar to the interpretation for the results in (3), the first and the second items are respectively the statistical error and the cost of privacy.

In the following, we present the minimax lower bound for the estimation in the high-dimensional Gaussian mixture model with differential privacy constraints, indicating the convergence rate obtained above is near optimal.

The proof of Proposition 4.3 is given in Section B.3. By comparing the results in Theorem 4.2 and Proposition 4.3, our algorithm is shown to attain the minimax optimality up to logarithm factors in the high-dimensional Gaussian mixture models.

Let $(\bm{x}_{1},y_{1}),(\bm{x}_{2},y_{2})...(\bm{x}_{n},y_{n})$ be the $n$ $i.i.d.$ observed samples from the mixture of regression model. Then, to use the EM algorithm, we need to compute

where $w_{\bm{\beta}}(\bm{x},y)=({1+\exp(-y\cdot\langle\bm{\beta},\bm{x}\rangle/{\sigma^{2}})})^{-1}.$

According to the gradient EM update rule, for the $t$-th iteration $\bm{\beta}^{t}$, we update $\hat{\bm{\beta}}$ by

Similar to the Gaussian Mixture model, to apply Algorithm 2, we need to specify the truncation operator $f_{T}(\nabla Q_{n}(\bm{\beta}^{t};\bm{\beta}^{t}))$. Rather than using truncation on the whole gradient, we perform the truncation on $y_{i}$, $\bm{x}_{i}$ and $\bm{x}_{i}^{\top}\bm{\beta}$ respectively, which leads to a more refined analysis and improved rate in the statistical analysis. Specifically, we define

Due to space constraints, we present the full algorithm for mixture of regression model in section A.1. We also verify Conditions 3.4-3.7 in the mixture of regression model, and summarize the results in Lemma A.1. In the following, we show the theoretical guarantees for the high-dimensional DP EM algorithm on the mixture of regression model, with proof given in Section B.4.