\subsection{Inference on Centroids}
Given the $K$ initial centroids $\vec{c}_1^{(0)}, \ldots,
\vec{c}_K^{(0)}$ (chosen independent of $D$),
let $\mathcal{S} = \langle\tilde{\theta}_1, \ldots, \tilde{\theta}_T\rangle$ be the
sequence of (noisy) outputs generated by running Algorithm~\ref{alg:dpkm} for $T$
iterations, where $\tilde{\theta}_t=(\tvec{s}_1^{(t)}, \ldots, \tvec{s}_{K}^{(t)},
  \tilde{n}_1^{(t)}, \ldots, \tilde{n}_K^{(t)})$ for $t=1,
\cdots, T$. Notice that (noisy) cluster centroids $\vec{c}_1^{(t)}, \ldots,
\vec{c}_K^{(t)}$ are completely determined by $\tilde{\theta}_t$. We
abuse notation and use $\tilde{\theta}_t$ to denote both noisy
statistics and $K$ centroids at iteration $t$.
Let $S(D,\theta)$ and $N(D,\theta)$ be the functions
that return the sum and the number of data points in each partition
determined by the given centroids $\theta$. 

Our goal is to make an inference on the cluster centroids
$\vec{c}_1, \ldots, \vec{c}_K$ based on the information $\mathcal{S}$
we learned privately from $D$. We do this by simulating datasets and
evaluating the likelihood of the observed noisy statistics
$\mathcal{S}$ under each dataset. 
Once a dataset that maximizes the likelihood of $\mathcal{S}$ is
found, new estimates for the cluster centroids can be derived by
running a non-private $K$-means algorithm (possibly with multiple random
restarts) on the dataset.
The log-likelihood
is defined by:
\begin{align*}
&\ln \mathbb{P}[\mathcal{S} \given D] 
=\ln \mathbb{P}[\tvec{s}_1^{(1)}, \ldots, \tvec{s}_K^{(T)}, \tilde{n}_1^{(1)}, \ldots,
  \tilde{n}_K^{(T)} \given D] \\
&=\sum\nolimits_{t=1}^T\left(\ln \pr\nolimits\p[\tvec{s}_1^{(t)},\ldots,\tvec{s}_K^{(t)} \given
  S(D,\tilde{\theta}_{t-1})]\right. \\
&\qquad\qquad + \left. \ln \pr\nolimits\p[\tilde{n}_1^{(t)}, \ldots,
  \tilde{n}_K^{(t)} \given   N(D,\tilde{\theta}_{t-1})]
  \right) \\
&\propto \sum_{t=1}^T\sum_{k=1}^K \norm*{\tvec{s}_k^{(t)} -
  S_k(D,\tilde{\theta}_{t-1})}_1 + \norm*{\tilde{n}_k^{(t)} -
  N_k(D,\tilde{\theta}_{t-1})}_1 \,,
\end{align*}
where the subscript $k$ in
$S_k(D,\tilde{\theta}_{t-1})$ and $N_k(D, \tilde{\theta}_{t-1})$
represent the sum and number of data points in the $k$\th cluster,
respectively.

\subsection{Imposing Consistency}
The accuracy of noisy output $\mathcal{S}$ can be improved by
imposing consistency constraints, using the algorithm proposed
in~\cite{Lee2015admm}. Suppose $\hvec{s}_k^{(t)}$ and
$\hat{n}_k^{(t)}$ are new estimates for $\tvec{s}_k^{(t)}$ and
$\tilde{n}_k^{(t)}$, respectively. It is clear that they should
satisfy the following constraints:
\begin{align*}
& \sum\nolimits_{k=1}^K \hvec{s}_k^{(1)} = \sum\nolimits_{k=1}^K \hvec{s}_k^{(2)} =
\ldots = \sum\nolimits_{k=1}^K \hvec{s}_k^{(T)}\,,\\
& \sum\nolimits_{k=1}^K \hat{n}_k^{(1)} = \sum\nolimits_{k=1}^K \hat{n}_k^{(2)} =
\ldots = \sum\nolimits_{k=1}^K \hat{n}_k^{(T)}\,,\text{ and} \\
& \hat{n}_k^{(t)} \geq 0 \text{ for all } k=1, \ldots, K \text{ and } t=1,
\ldots, T\,.
\end{align*}
For clear semantics and better readability, in the following we
continue to use the notation $\tvec{s}_k^{(t)}$ 
and $\tilde{n}_k^{(t)}$, but they represent the post-processed values.

\subsection{Simulation via MCMC}
The proposed algorithm makes use of approximate sampling method to
find a dataset under which the likelihood of $\mathcal{S}$ is 
maximized. Using MCMC, it samples datasets from the approximate
posterior distribution $\pr\nolimits\p(D \given \mathcal{S})$ and evaluates the
likelihood, while keeping track of the best solution. The target
distribution is 
\[
\begin{aligned}
\pi(D) &= \Exponential\left(-\sum\nolimits_{t=1}^T\sum\nolimits_{k=1}^K
\epsilon_{\!_{S}}\norm{\tvec{s}_k^{(t)} - 
  S_k(D,\tilde{\theta}_{t-1})}_1\right.\\
&\qquad \qquad + \left.\epsilon_{\!_{N}}\norm{\tilde{n}_k^{(t)} -
  N_k(D,\tilde{\theta}_{t-1})}_1\right)\,,
\end{aligned}
\]
where $\epsilon_{\!_{S}}$ and $\epsilon_{\!_{N}}$ correspond to the
privacy budgets for noisy cluster sums and sizes.\footnote{For
  simplicity, we assume $\epsilon_{\!_{S}} = \epsilon_{\!_{N}}$.}
\paragraph*{Proposal distribution}
The Metropolis-Hastings (MH) algorithm requires choice of proposal
distribution, and the convergence of Markov chain to its stationary
distribution $\pi$ is greatly dependent on that choice. The use of a
proposal distribution that is far from $\pi$
will have a high rejection rate and result in slow convergence.

It is shown that $K$-means algorithm can be thought as a special case of
Gaussian Mixture Model (GMM), with means equal to centroids and a common
covariance set to $\delta\mat{I}$ for small $\delta>0$. Given
$\theta_t = (\tvec{s}_1^{(t)}, \ldots,
\tvec{s}_K^{(t)},\tilde{n}_1^{(t)}, \ldots, \tilde{n}_K^{(t)})$, our
proposal distribution is defined to be a mixture of Gaussians:
\begin{equation}
q(\vec{x}) = \sum_{k=1}^K \omega_k\Norm{\vec{x}; \vec{c}_{k}^{(t)},
  \delta\mat{I}}\,,
\label{eq:proposal}
\end{equation}
where $\omega_k = \slfrac{\tilde{n}_k^{(t)}}{\sum_{i=1}^K
  \tilde{n}_i^{(t)}}$ and $\vec{c}_k^{(t)} = \tvec{s}_{k}^{(t)} /
\tilde{n}_k^{(t)}$.  

Given the current dataset $D^{(\tau)}$, a new dataset $D'$ is proposed
by randomly choosing a data point $\vec{x}_i$ from $D$ and replacing
it with a new point $\vec{x}'$ sampled from the proposal
distribution $q$. The sampling of $\vec{x}'$ is done as follows:
\begin{enumerate}[label=(\roman*),leftmargin=*,topsep=1pt,itemsep=-.5ex]
\item choose an integer $t$ randomly from $\{1, 2, \ldots, T\}$.
\item sample $z\given t \sim \Cat{K, \omega_1, \ldots,
    \omega_K}$.
\item sample $\vec{x}'\given z, t \sim \prod_{k=1}^K
  \Norm{\vec{c}_k^{(t)}, \delta\mat{I}}^{\I{}{z=k}}$. 
\end{enumerate}
In the above, $\Cat{K, \omega_1, \ldots, \omega_K}$ represents the
Categorical distribution having possible values in $\{1, \ldots, K\}$,
each with probability mass $\omega_k$ for $k=1,\ldots, K$. $\I{}{z=k}$
is an indicator function whose value is 1 if $z=k$ and 0 otherwise.
\paragraph*{MH Algorithm}
Given $D^{(\tau)}$, the proposed dataset $D'$ of next state
$(\tau+1)$ is accepted with probability
\[
A(D^{(\tau)}, D') = \min\left\{\frac{\pi(\mathcal{S}\given
    D')}{\pi(\mathcal{S}\given D^{(\tau)})}
    \frac{q(\vec{x} \given \vec{x}')}{q(\vec{x}' \given \vec{x})},
    1\right\}\,.
\]
Without loss of generality, suppose a data point $\vec{x}$ is removed
from the $i$\th cluster and a new data point $\vec{x}'$ is added to
the $j$\th cluster at time $\tau$. Then we have
\begin{align*}
&\ln\pi(\mathcal{S}\given D') -\ln\pi(\mathcal{S} \given D^{(\tau)}) \\
&= -\sum_{t=1}^T
\norm{\tvec{s}_i^{(t)} - S_i(D',\tilde{\theta}_{t-1})}_1 -
\norm{\tvec{s}_i^{(t)} - S_i(D^{(\tau)},\tilde{\theta}_{t-1})}_1 \\
& \quad +
\norm{\tvec{s}_j^{(t)} - S_j(D',\tilde{\theta}_{t-1})}_1 - 
\norm{\tvec{s}_j^{(t)} - S_j(D^{(\tau)},\tilde{\theta}_{t-1})}_1 \\
& \quad + 
\abs{\tilde{n}_i^{(t)} - N_i(D',\tilde{\theta}_{t-1})} -
\abs{\tilde{n}_i^{(t)} - N_i(D^{(\tau)},\tilde{\theta}_{t-1})} \\
& \quad + 
\abs{\tilde{n}_j^{(t)} - N_j(D',\tilde{\theta}_{t-1})} -
\abs{\tilde{n}_j^{(t)} - N_j(D^{(\tau)},\tilde{\theta}_{t-1})}\,. 
\end{align*}
We note that $S(D',\tilde{\theta}_{t-1})$ and $N(D',
\tilde{\theta}_{t-1})$ can be calculated from $S(D^{(\tau)},
\tilde{\theta}_{t-1})$ and $N(D^{(\tau)}, \tilde{\theta}_{t-1})$,
respectively. 
The MH correction term is given by
\begin{align*}
\frac{q(\vec{x})}{q(\vec{x}')} 
&= \frac{\qq(\vec{x} \given \vec{z},
  t)\qq(\vec{z}=i)\qq(t)}{\qq(\vec{x}'\given \vec{z}',
  t')\qq(\vec{z}'=j)\qq(t')}  
= \frac{\tilde{n}_i^{(t)}\mathcal{N}(\vec{x}; \vec{c}_i^{(t)},
    \delta\mat{I})}{\tilde{n}_j^{(t')}\mathcal{N}(\vec{x}';\vec{c}_j^{(t')}, \delta\mat{I})}\,.
\end{align*}
The initial state $D^{(0)}$ is initialized with $K$ centroids at the
last iteration. It consists of $\tilde{n}_k^{(T)}$ data points at
$\vec{c}_{k}^{(T)}$ for $k=1, \ldots, K$.





\section{Introduction}
\label{sec:intro}
\input{intro}

\section{Related Works}
\label{sec:related}
\input{related}

\section{Postprocessing for K-means}
\label{sec:algorithm}
\input{algorithm}

\section{Experiments}
\label{sec:experiment}
\input{experiment}

\section*{Acknowledgement}
This research was supported by NSF grant 1228669.



