
\subsection{Datasets}
We used 4  large-scale datasets for evaluation.

\textbf{\underline{Partially synthetic housing}}. Individuals live in households and group quarters. The number of individuals in each facility is important but this information was truncated past households of size 7 in the 2010 Decennial Census, Summary File 1 \cite{sf1}. Thus we created a partially synthetic dataset that mirrors the published statistics, but adds a heavy tail as would be expected from group quarters (e.g., dormitories, barracks, correctional facilities). This was done for each state by estimating the ratio $\#$ households of size 7/$\#$ households of size 6, and then randomly sampling (with a binomial distribution) the number of groups of size $k\geq 8$ so that the same ratio holds (in expectation) between number of groups of neighboring sizes. Then 50 outliers groups are chosen with size uniformly distributed between 1 and 10,000. The hierarchy in this levels are National and State (50 states plus Puerto Rico and the District of Columbia). The third level is County, which we obtained by randomly assigning groups at the state to their counties (the assignment was proportional to county size).

\textbf{\underline{NYC taxi}}: We use 143,540,889 Manhattan taxi trips from the 2013 New York City taxi dataset \cite{nyctaxi}. An anonymized taxi medallion (e.g., a taxi) is considered a group and the size of the group is the number of pickups it had in a region. The region hierarchy is the following. Level 0: Manhattan; level 1: upper/lower Manhattan; level 2: 28 neighborhoods from NTA boundary \cite{nycopen}.

\textbf{\underline{Race distribution}} (white and Hawaiian): For each block, based on 2010 Census data (in Summary File 1 \cite{sf1}), we count the number of whites and number of native Hawaiians that live in the block. Hence block is treated as a group. The hierarchy is National, State, and County. We performed evaluations on all 6 major race categories recorded by the Census, but omitted the rest due to space restrictions. 

The statistics for our datasets are the following.\\
\begin{tabular}{ | c | c | c | c |  }
\hline
   Data      &   \# groups  & \# people/trip  &   \# unique size  \\
\hline 
Synthetic & 240,908,081 & 605,304,918 & 2352 \\
\hline
White & 11,155,486 & 226,378,365 & 1916 \\
\hline
Hawaiian & 11,155,486 & 540,383 & 224\\
\hline
Taxi & 360,872 & 130,962,398 & 3128 \\
\hline
\end{tabular}

\begin{comment}
We evaluate our algorithm on 3 types of datasets with two hierarchies. 
  
 The datasets include one partially synthetic household size dataset and the rest are real-world data: two race data, taxi trip data. 
  
The two hierarchies are a \textbf{2-level hierarchy} with level 0 as the United States, level 1 as 52 states in the entire country and a \textbf{3-level hierarchy} with level 0 as the west coast of the United States, level 1 as 3 states of the west coast - California, Oregon, Washington state, level 2 as 133 counties which belong to the states at west coast. In the context of taxi trips, the 3-level hierarchy represents level 0 as manhattan, level 1 as the upper/lower town, level 2 as the neighborhoods.
\end{comment}
For some \changes{count-of-counts\xspace} estimation methods, such as the cumulative sum $H_c$ method, one needs to specify a public maximum group size $K$. We set $K=100,000$ as a conservative estimate (for example, in our partially synthetic housing dataset, the true max size was around $10,000$, an order of magnitude smaller). 
For our 2-level hierarchy experiments on Census related data, we used National/State. For 3 level, we used West Coast/State/County. All numbers plotted are averaged over 10 runs.





\begin{figure*}[ht]
	\centering
	\includegraphics[width=14.5cm]{myfig1_reduce_space.png}
 \vspace{-7pt}
	\caption{Merging estimates using weighted average vs. normal average. x-axis: privacy budget per level.}
  \label{fig:metrics}
\end{figure*}



\subsection{Evaluation}
\textbf{Evaluation metric}. For each level of the hierarchy, we evaluate Earthmover's distance (emd) as discussed in \Cref{subsec:error}, per node in the level, in order to see error at each level. We do not aggregate error across all levels of the hierarchy as there is no principled way of, for example, weighting the importance of error at the state level compared to the national level. Each error measure is averaged across 10 runs. The standard deviation of the average is then the empirical standard deviation (which gives std for one random run) divided by $\sqrt{10}$ (which then gives the std for the mean of 10 runs). We plot 1 std error bars on each figure.

\textbf{Algorithm selection}. We evaluate a variety of choices for generating hierarchical \changes{count-of-counts\xspace histograms}. For example, we can use the $H_g$ method to generate estimates at national level but $H_c$ at state level and $H_c$ again at county level, which we would denote as $H_g\times H_c\times H_c$ and compare it against $H_g\times H_g\times H_g$. \changes{In general, using the $H_c$ method at each level gives best results except on the partially synthetic data, thus we recommend it as the default choice. However, as is true with all data-dependent, no algorithm is expected to dominate on all datasets \cite{Kotsogiannis2017:Pythia}. For fine-grained algorithm selection, one could use Pythia \cite{Kotsogiannis2017:Pythia} or the approach of Chaudhuri et al. \cite{Chaudhuri2011ERMJMLR}. As these are well-established techniques, they are outside the scope of the paper.}

\textbf{Interpreting error}. How can we interpret the error numbers of the algorithms? I.e., what is a good error? We suggest comparison against the following ``omniscient'' algorithm. Given an $\epsilon$, the omniscient algorithm will know which group sizes exist in which node in the hierarchy and reduce it to the problem of estimating a simple histogram of known group size by location. Thus it will split its privacy budget per level of the hierarchy and will add Laplace$(1/\epsilon)$ noise (with standard deviation of $\sqrt{2}/\epsilon$) only to those groups that exist. Meanwhile, the $\epsilon$-differentially private algorithms must effectively estimate which group sizes exist and also estimate their sizes, along with which nodes in the hierarchy the groups belong to. The error for the omniscient algorithm is then expected to be $[\#\text{distinct group sizes} \times \sqrt{2}/\epsilon\times \#\text{levels}]$. For example, in Figure \ref{fig:metrics}, at privacy budget $0.1$ per level and the partially synthetic housing data at the national level, there are $2,352$ groups and so the omniscient algorithm will have expected error around $3.3\times 10^4$, which is in line with the results for the best differentially private method in that figure.









\subsubsection{Ruling out the naive method}\label{subsec:expnaive}
The naive strategy adds noise to $\tau.H$ directly (Section \ref{subsec:naive}).   Its average error with $\epsilon=1$ is \changes{in the billions, as} shown in the following table \changes{confirming the analysis in \Cref{subsec:naive}.} 

\smallskip


\hspace{-14pt}
\begin{tabular}{ | c | c | c | c |  }
\hline
 Synthetic    &  White  & Hawaiian  & Taxi  \\
\hline
4,462,728,374 &  4,809,679,734 & 4,027,891,692 & 208,977,518\\
\hline
\end{tabular}
\changes{This error is several orders of magnitude larger than the other methods (e.g. Figure \ref{fig:metrics}), so is not considered further.}





\begin{figure*}[htb!]
	\centering
	   
    \includegraphics[width=14.5cm]{myfig2_reduce_space.png}
\vspace{-5pt}
	\caption{2-level consistency at each level. x-axis: privacy budget per level.}
  \label{fig:2lv-all}
\end{figure*}


\subsubsection{Comparison to bottom-up aggregation}
\changes{
The Bottom-Up (BU) baseline allocates all privacy budget to the leaves and then sets the count-of-counts\xspace histogram of a parent to be the sum of the histograms at the leaves. As mentioned earlier, it is expected to have very low error at the leaves but higher error everywhere else (an observation that is true in general for hierarchical problems \cite{Qardaji2013Hierarchical,Hay2010Boosting}). In the following table we examine the error at each level of BU compared to the consistency algorithm that uses $H_c$ at each level. Both approaches use a total privacy budget of $\epsilon=1.0$.
}

\vspace{-0.2em}
\begin{tabular}{c|r|r|r|r|}\cline{2-5}
\multicolumn{1}{c}{} & \multicolumn{1}{|c}{Part. Synth.} & \multicolumn{1}{|c}{White}& \multicolumn{1}{|c}{Hawaiian} & \multicolumn{1}{|c|}{Taxi}\\\cline{2-5}
\multicolumn{0}{c}{}&\multicolumn{4}{c}{Level 0}\\\cline{2-5}
BU &  $78,459.0$ & $448,909.0$ & $13,968.0$ & $20,731.0$\\
$H_c$ & $32,480.0$ & $17,000.0$ & $1,381.0$ & $10,547.0$\\\cline{2-5}
\multicolumn{0}{c}{}&\multicolumn{4}{c}{Level 1}\\\cline{2-5}
BU &  $1,512.2$ & $8,722.0$ & $270.1$ & $10,405.5$\\
$H_c$ & $ 1,000.3$ & $1,511.8 $ & $ 117.7$ & $ 5,431.5$\\\cline{2-5}
\multicolumn{0}{c}{}&\multicolumn{4}{c}{Level 2}\\\cline{2-5}
BU &  $24.9 $ & $152.3$ & $4.3$ & $772.8$\\
$H_c$ & $ 80.1$ & $363.8$ & $21.6$ & $1,601.8$\\\cline{2-5}
\end{tabular}

\changes{The results are as expected, with significant improvements at the higher level nodes in exchange for a small increase in error at the leaves.}


\subsubsection{Weighted Average Estimation Comparison} \label{subsec:wavg}

Next we evaluate the merging strategy for groups that are matched across different levels of the hierarchy (as described in Section \ref{subsec:merge}). There are two choices. If a group from one level of the hierarchy (e.g., a household at the state level) is matched to a group at a different level (e.g., at the national level), we can either average the two group size estimates, or use a weighted average based on their estimated variances from Section \ref{subsec:initvar}. If the variance estimates are accurate, we expect such a merging to result in improved histogram error when compared to normal averaging (see Figure \ref{fig:metrics}).



We estimate two levels of the hierarchy for each dataset (partially synthetic, Hawaiian race distribution, White race distribution are shown due to space limitations) and consider various methods for the initial group size estimates in each of the two levels (Section \ref{sec:nonhier}); for example $H_c\times H_g$ means that the $H_c$ method was used for the top level and $H_g$ for the second level. We see that the weighted average method consistently produces large reductions in error (compared to normal averaging) at the top level for every privacy budget and combination of methods ($H_g/H_c$), and consistently produces modest improvements in error at the second level. Hence we conclude that our variance estimates, which are used in the merging procedure, are useful approximations. The only method missing from these graphs is the combination $H_g\times H_g$. On these datasets, the error of normal average was so large that it would visually skew the results. Due to the superiority of weighted averaging, all following experiments use this method of merging estimates.



\begin{comment}
\begin{figure}[ht!]
	\centering
	\includegraphics[width=8cm]{figs/Hcwhite-s0s2000.png}
	\caption{$x$-axis: group sizes. $y$-axis: the group count in log scale. Top: estimated size histogram. Bottom: true size histogram}
  \label{fig:dist}
\end{figure}
\end{comment}







\subsubsection{2-Level Hierarchy Results} \label{subsec:2lv}


We next show our results for the case when  2 levels of the hierarchy need to be estimated, in  \Cref{fig:2lv-all}. Due to space limitations, we present the consistency result for $H_g\times H_g$ and $H_c\times H_c$ using weighted averaging to merge estimates. One thing to notice is that the best performing method is comparable to the omniscient baseline \changes{and generally the $H_c$ method is the one that performs best.} Typically, we expect the $H_c$ method to dominate for ``dense'' data. For example, white population data contains many groups from size 0 to size 3000. On the other hand, the partially synthetic housing data is more sparse at the national level, with many small groups (e.g., sizes 1-12) followed by large gaps between group sizes. In such cases, methods based on $H_g$ work better.




 














\begin{comment}
\begin{figure*}[ht]
\centering
\subfigure[synthetic]{\includegraphics[width=0.3\textwidth]{figs/synth-woc.png}}
\subfigure[asian]{\includegraphics[width=0.3\textwidth]{figs/asian-woc.png}}
\subfigure[black]{\includegraphics[width=0.3\textwidth]{figs/black-woc.png}}
\subfigure[hawaiian]{\includegraphics[width=0.3\textwidth]{figs/haw-woc.png}}
\subfigure[native]{\includegraphics[width=0.3\textwidth]{figs/native-woc.png}}
\subfigure[white]{\includegraphics[width=0.3\textwidth]{figs/white-woc.png}}
\caption{2-level without consistency}
\label{fig:woc}
\end{figure*}



\begin{figure*}[ht]
\centering
\subfigure[synthetic]{\includegraphics[width=0.3\textwidth]{figs/synth-wc.png}}
\subfigure[asian]{\includegraphics[width=0.3\textwidth]{figs/asian-wc.png}}
\subfigure[black]{\includegraphics[width=0.3\textwidth]{figs/black-wc.png}}
\subfigure[hawaiian]{\includegraphics[width=0.3\textwidth]{figs/haw-wc.png}}
\subfigure[native]{\includegraphics[width=0.3\textwidth]{figs/native-wc.png}}
\subfigure[white]{\includegraphics[width=0.3\textwidth]{figs/white-wc.png}}
\caption{2-level after consistency}
\label{fig:wc}
\end{figure*


\begin{figure*}[ht]
\centering
\subfigure[synthetic-national]{\includegraphics[width=0.22\textwidth]{figs/synth-natall.png}}
\subfigure[white-national]{\includegraphics[width=0.22\textwidth]{figs/white-natall.png}}
\subfigure[hawaiian-national]{\includegraphics[width=0.22\textwidth]{figs/haw-natall.png}}
\subfigure[black-national]{\includegraphics[width=0.22\textwidth]{figs/black-natall.png}}
\subfigure[synthetic-state]{\includegraphics[width=0.22\textwidth]{figs/synth-stateall.png}}
\subfigure[white-state]{\includegraphics[width=0.22\textwidth]{figs/white-stateall.png}}
\subfigure[hawaiian-state]{\includegraphics[width=0.22\textwidth]{figs/haw-stateall.png}}
\subfigure[black-state]{\includegraphics[width=0.22\textwidth]{figs/black-stateall.png}}
\caption{2-level w. consistency v.s w.o. consistency at each level}
\label{fig:level}
\end{figure*}

\end{comment}


\subsubsection{3-Level Hierarchy Results}  \label{subsec:3lv}




\begin{figure*}[ht]
	\centering
	\includegraphics[width=18cm]{myfig3_reduce_space.png}
	\caption{3-level consistency at each level. x-axis: privacy budget per level.}
  \label{fig:3lv-all}
\end{figure*}



We next show the results when three levels of the hierarchy need to be estimated. Because there are over 3,000 counties (hence 3,000 isotonic regressions), for computational reasons we limit the hierarchy to the west coast (for the partially synthetic housing data as well as the race distribution data.). Of these two alternatives $H_g \times H_g \times H_g$ and $H_c \times H_c \times H_c$ on synthetic, race and taxi at each level in \Cref{fig:3lv-all}. The taxi data uses its full geography. In general, we see that no method dominates the other, but group size estimation methods based on $H_c$ generally perform better than $H_g$ and so are  a good default choice. 












\begin{comment}
\begin{figure*}[ht]
\centering
\subfigure[synthetic-West Coast]{\includegraphics[width=0.3\textwidth]{figs/synth-wcall.png}}
\subfigure[synthetic-West Coast States]{\includegraphics[width=0.3\textwidth]{figs/synth-wcstateall.png}}
\subfigure[synthetic-West Coast Counties]{\includegraphics[width=0.3\textwidth]{figs/synth-wcctall.png}}
\subfigure[white-West Coast]{\includegraphics[width=0.3\textwidth]{figs/white-wcall.png}}
\subfigure[white-West Coast States]{\includegraphics[width=0.3\textwidth]{figs/white-wcstateall.png}}
\subfigure[white-West Coast Counties]{\includegraphics[width=0.3\textwidth]{figs/white-wcctall.png}}
\caption{3-level w. consistency v.s w.o. consistency at each level}
\label{fig:level}
\end{figure*}

\end{comment}
  


\section{Introduction}\label{sec:intro}
\input{intro}
\section{Related Work}\label{sec:related}
\input{related}
\section{Problem Definition}\label{sec:notation}
\input{notation}

\section{The Base Case: Non-hierarchical Count-of-Counts}\label{sec:nonhier}
\input{nonhier}
\section{Hierarchical Consistency}\label{sec:consistency}
\input{method}
\section{Experiments}\label{sec:experiments}
\input{experiments}
\section{Conclusions and Future Work}\label{sec:future}
In this paper we introduced the differentially private \emph{hierarchical} \changes{count-of-counts\xspace} histogram problem and presented a solution based on isotonic regression and optimal weighted matchings. This problem is motivated by a variety of tables that are published in truncated form in Summary File 1 of the 2010 U.S. Decennial Census. The actual tables include additional demographic characteristics that are attached to the household sizes at each level of geography. Such additional information greatly expands the dimensionality of the problem and is an area of future work.



\smallskip
\section{Acknowledgments}
We thank John Abowd for discussions on disclosure avoidance requirements for publishing Census tables. \changes{We thank the anonymous reviewers for their help comments}. The work was partially supported by NSF grants 1054389, 1544455, and 1702760. The views in this paper are those of the authors, and do not necessarily represent the Census Bureau.
\balance

\clearpage
\bibliographystyle{abbrv}

\subsection{Initial Variance Estimation}\label{subsec:initvar}
The first part of our algorithm produces a differentially private \changes{count-of-counts\xspace} histogram $\tau.\widehat{H}$ for every node $\tau$. We then convert it into the \changes{unattributed} histogram $\tau.\widehat{H}_g$. For each  $i$, we need an estimate of the variance of the $i^\text{th}$ largest group $\tau.\widehat{H}_g[i]$. Now, the \changes{count-of-counts\xspace} histogram $\tau.\widehat{H}$ is obtained either from the  $H_g$ method (Section \ref{subsec:unatt}) or the cumulative histogram $H_c$ method (Section \ref{subsec:cumu}) and so the variance depends on which method was used.\footnote{\changes{Generally, $H_c$ works well for all levels. Users preferring fine-grained control can use generic algorithm selection tools \cite{Kotsogiannis2017:Pythia,Chaudhuri2011ERMJMLR}.}}

\subsubsection{Variance estimation for the $H_g$ method} \label{subsubsec:HgVar}
Recall that in the \changes{unattributed} histogram method, we obtained a noisy array (noise was added to $\tau.H_g$) and performed isotonic regression on it to get $\tau.\widehat{H}_g$. Since we need an estimate of the variance  of $\tau.\widehat{H}_g[i]$, we can use the following special properties of isotonic regression solutions  \cite{barlowstatistical}. As shown in Figure \ref{fig:l2}, isotonic regression is equivalent to taking a noisy 1-d array, partitioning the array, and assigning the same value to each element within a partition. In the case of $L_2$ isotonic regression, this value is the average of the noisy counts from the partition. In the case of $L_1$, this value is the median of the noisy counts in the partition. Thus the variance of each cell in the resulting array depends on the variance due to partitioning and the variance due to averaging the noisy counts. We cannot quantify the variance due to partitioning. However, the variance due to averaging noisy counts can be estimated. 

Each noisy count is generated by adding noise from the double-geometric distribution with scale $1/\epsilon_1$. Its variance can be approximated by the variance of the Laplace distribution with the same scale. Namely, the variance is $2/\epsilon_1^2$. In a partition of size $S$ (and $L_2$ isotonic regression), the value assigned to that partition is the average\footnote{In the case of $L_1$ regression, the value assigned to the partition is the median of $S$ noisy values; the variance of the median is difficult to compute, so we again estimate it as $2/(S\epsilon^2)$} of $S$ noisy values and so has variance $2/(S\epsilon^2)$. The partitions that were created by the solution of the isotonic regression are easy to determine -- they are simply the consecutive entries in the solution that have the same value.

Thus, our estimate of the variance for group $i$ in $\tau.\widehat{H}_g$ is computed as follows. Let $S_i$ be the number of groups that were in the same partition as $i$ in the solution (i.e. the number of entries in $\tau.\widehat{H}_g$ that equal $\tau.\widehat{H}_g[i]$). Set the variance estimate for the $i^\text{th}$ largest group to be  the following in Line \ref{line:estvar} of Algorithm \ref{alg:top-down}: $\tau.V_g[i]= \frac{2}{|S_i|\epsilon_1^2}.$

\begin{figure}[t!]
\centering
\includegraphics[width=0.28\textwidth]{isotonic.pdf}
\vspace{-9pt}
\caption{$L_2$ Isotonic Regression converts a noisy histogram that is no longer nondecreasing into a nondecreasing histogram by partitioning and averaging within each partition.}
\label{fig:l2}
\end{figure}




\begin{comment}
\begin{figure}
\begin{align*}
\toprule
~~~~~\tau.H = ~ [~0, ~0, &~1, ~2, ~3~] \\
~~~~~\tau.H_g = ~ [~2, ~3, &~3, ~4, ~4, ~4~] \\
~~~~~\tau.\widetilde{H_g} =   [~0, ~4, &~2, ~4, ~5, ~3~] \\
				   &	\Downarrow L_2  \text{~isotonic regression} \\
~~~~~\tau.\widetilde{H_g} =   [~0, ~~||~~ \underbrace{4, ~2}, & ~~||~~\underbrace{4,~ 5, ~3}~]  \textbf{~~partition~~}\\	
~~~~~\tau.\widehat{H_g} =   [~\underbrace{0,} ~~||~~ \underbrace{3, ~3}, & ~~||~~\underbrace{4,~ 4, ~4}~]  \textbf{~~taking average~~}\\
~~~~~~~~~~~~~~~~~~~~~~~~~~S_0	~~~~~~~~S_1	~~ & ~~~~~~~~~  S_2 ~~~~~~~~~~~\\		 							 
\bottomrule
\end{align*}
\vspace{-30pt} 
\caption{$L_2$ Isotonic Regression}
\label{fig:l2}
\end{figure}
\end{comment}


\subsubsection{Variance for the cumulative histogram method} \label{subsubsec:HcVar}
The cumulative histogram method also uses isotonic regression, but uses it to create $\tau.\widehat{H}_c$ while we need estimates of the variance of $\tau.\widehat{H}_g[i]$. Since the $H_g$ representation is a nonlinear transformation of the $H_c$ representation, we must use a different way to estimate variance.

 We know that before the isotonic regression in the $H_c$ method, independent noise with scale $1/\epsilon_1$ was added to each cell of $\tau.H_c$ so we (over)estimate the variance of $\tau.\widehat{H}_c[j]$ as $2/\epsilon_1^2$. The estimated number of groups of size $j$ is $\tau.\widehat{H}_{c}[j]-\tau.\widehat{H}_{c}[j-1]$ so it would have variance $4/\epsilon_1^2$. Dividing by the number of estimated groups of size $j$, this gives $4\epsilon_1^{-2}/(\tau.\widehat{H}_{c}[j]$ $-\tau.\widehat{H}_{c}[j-1])$ as the estimated variability in each group whose size is estimated to be $j$. Thus we set our estimate of the variance of the $i^\text{th}$ largest group to be the following in Line \ref{line:estvar} of Algorithm \ref{alg:top-down}:
$$\textstyle{\tau.V_g[i]=4 / (\epsilon_1^2\times\text{number of estimated groups of size $\tau.\widehat{H}_g[i]$})}$$


\subsection{Optimal Matching}
\label{subsec:match}
Now, every group belongs to a node in each level of the hierarchy (e.g., a household in Fairfax county is also a household in Virginia). Since group size distributions are estimated independently at all levels of the hierarchy (lines \ref{line:tdlevel} -- \ref{line:tdswitch} in Algorithm \ref{alg:top-down}), this causes several problems. First, each group has several size estimates (a household has a certain estimated size from the Fairfax County estimates and another estimated size from the Virginia estimate). Second, from these separate estimates, we don't know which group in one level of the hierarchy corresponds to which group at a different level of the hierarchy (e.g., is it possible that the $102^\text{nd}$ largest household in Virginia, with estimated size 12 is the same as the $3^\text{rd}$ largest household in Fairfax County, with estimated size 13?).

Thus, to make the group size \changes{distribution} estimates consistent, we need to estimate a matching between groups in $\tau$ and groups in the children of $\tau$, as in Figure \ref{fig:illustration}.
For privacy reasons, such a matching must only be done using the differentially private data generated in Line \ref{line:tdswitch} of Algorithm \ref{alg:top-down}. In this section we explain how to perform this matching and in the next section we explain how to reconcile the different size estimates.

Formally, a matching is a function $m$ that inputs a node $\tau$ and an index $i$ and returns a node $c$ and index $j$ with the semantics that the $i^\text{th}$ smallest group in $\tau$ is believed to be the same as the $j^\text{th}$ smallest group in child $c$ of $\tau$. We must estimate this matching using only differentially private data (e.g., $\tau.\widehat{H}$ for each node $\tau$). We first convert each $\tau.\widehat{H}$ into the \changes{unattributed} representation $\tau.\widehat{H}_g$, where $\tau.\widehat{H}_g[i]$ is the size of the $i^\text{th}$ smallest group in $\tau$.

For each node $\tau$, we set up a bipartite weighted graph as shown in Figure \ref{fig:illustration}. There are $\tau.G$ nodes on the top half. We label them as $(\tau, 1), (\tau,2),\dots, (\tau, \tau.G)$. Each node on the bottom half has the form $(c, j)$, where $c$ is a child of $\tau$ and $j$ is an index into $c.\widehat{H}_g$. Between every node $(\tau, i)$ and $(c, j)$ there is an edge with weight $|\tau.\widehat{H}_g[i]-c.\widehat{H}_g[j]|$ which measures the difference in estimated size between the $i^\text{th}$ smallest group in $\tau$ and the $j^\text{th}$ smallest group in $c$. 

Our desired matching is then the least cost weighted matching on this bipartite graph. Sophisticated matching algorithms (e.g., based on network flows) can find an optimal matching, but they have time complexity at least $O(\tau.G)^3$ \cite{lovasz}. In our case, $\tau.G$ can be in the millions (e.g., there are over 100 million households in the U.S.).

There exists a well-known 2-approximation algorithm for matching, which adds edges to the matching in order of increasing weight \cite{lovasz}. However on our graph it would run in $O(\tau.G^2\log\tau.G)$ time as there are $\tau.G^2$ edges and they need to be sorted. Instead, we take advantage of the special weight structure of our edges to produce an optimal matching algorithm with time complexity $O(\tau.G\log\tau.G)$.

\begin{figure}[t!]
\centering
\includegraphics[width=0.38\textwidth]{matching-new.pdf}
\vspace{-9pt}
\caption{Consistency Matching Illustration}
\label{fig:illustration}
\end{figure}



\subsubsection{Optimal Matching Algorithm}
\input{newalg}

\subsection{Merging Estimates}\label{subsec:merge}
Given a node $\tau$, the matching algorithm assigns one group in $\tau$ to one group in some child of $\tau$ (i.e. it says that, for every $i$, the $i^\text{th}$ smallest group in $\tau$ matches the $j^\text{th}$ smallest group in child $c$ of $\tau$). This means that for every group, we have two estimates of its size: $\tau.\widehat{H}_g[i]$ and $c.\widehat{H}_g[j]$ as well as corresponding estimates of its variance $\tau.V_g[i]$ and $c.V_g[j]$.
There are two possible ways of reconciling these estimates.

\textbf{Naive strategy.}
The simplest way is to simply average  $\tau.\widehat{H}_g[i]$ and $c.\widehat{H}_g[j]$. This approach would be valid if the variance estimates were not accurate -- recall that it is not possible to estimate the variances exactly, so they needed to be approximated. However, as we show in our experiments, a weighted averaging based on the estimated variance outperforms this strategy.

\textbf{Variance-based weighted strategy.} If we have two noisy estimates of the same quantity (e.g., $\tau.\widehat{H}_g[i]$ and $c.\widehat{H}_g[j]$), along with their variances  $\tau.\widehat{H}_g[i]$ and $c.\widehat{H}_g[j]$, it is a well-known statistical fact (e.g., see \cite{Hay2010Boosting}) that the optimal linear way of combining the estimates is to estimate the size as a weighted average, where the weights are inversely proportional to the variance:
{
\begin{align}
\textstyle{
\left(\frac{\tau.\widehat{H}_g[i]}{\tau.V_g[i]} + \frac{c.\widehat{H}_g[j]}{c.V_g[j]}\right)  \Big/\left(\frac{1}{\tau.V_g[i]} + \frac{1}{c.V_g[j]}\right)\label{eqn:varest1}
}
\end{align}
}

and the variance of this estimator is
{
\begin{align}
\textstyle{
\left(\frac{1}{\tau.V_g[i]} + \frac{1}{c.V_g[j]}\right)^{-1}\label{eqn:varest2}
}
\end{align}
}
Thus we would update the size of the $j^\text{th}$ largest group at child $c$ using Equation \ref{eqn:varest1} and its new variance as Equation \ref{eqn:varest2}. This estimator is preferable to the naive strategy when the variance estimates are accurate.

The size estimates are then rounded. After the size estimates at the children are updated, the top-down algorithm continues matching the groups in each child $c$ with the groups at the children of $c$. Once the groups at the leaves are updated, the resulting sizes at the leaves are treated as the final estimates.

 
\subsection{Privacy}
\begin{theorem}
Algorithm \ref{alg:top-down} satisfies $\epsilon$-differential privacy.
\end{theorem}
\begin{proof}
Privacy is easy to analyze because the algorithm separates the differentially private data access from the post-processing. Specifically, the only part of the algorithm that touches the sensitive data occurs in Lines \ref{line:tdlevel} through \ref{line:tdswitch}. It uses sequential composition across levels of the hierarchy. Thus each of the $L+1$ levels is assigned $\epsilon/(L+1)$ of the privacy budget. Within each level there is parallel composition because adding or removing one person from a group only affects the node that contains that group (and none of the sibling nodes). The \changes{count-of-counts\xspace} histograms produced at each node use either the method of Section \ref{subsec:unatt} or Section \ref{subsec:cumu}) with privacy budget $\epsilon/(L+1)$ and they scale the noise correctly to the global sensitivity, as discussed in those sections.

The rest of the top-down algorithm is completely based on the differentially private results of of Lines \ref{line:tdlevel} through \ref{line:tdswitch}. The conversions between $\tau.\widehat{H}_g, \tau.\widehat{H}_c,$ and $\tau.\widehat{H}$ are trivial manipulations of the histogram format (they do not touch the original data). The variance estimation is based on $\tau.\widehat{H}_g, \tau.\widehat{H}_c,$ and $\tau.\widehat{H}$. The matching algorithm only uses $\tau.\widehat{H}_g$ for each node, and the method for merging estimates only uses $\tau.\widehat{H}_g$ (for each node) and the associated variance estimates (which were computed from $\tau.\widehat{H}_g, \tau.\widehat{H}_c,$ and $\tau.\widehat{H}$). Since post-processing differentially private outputs still satisfies differential privacy \cite{DMNS06Calibrating}, the overall algorithm satisfies differential privacy.
\end{proof}










 



\subsection{Naive Strategy}\label{subsec:naive}
The naive strategy for estimating $\widehat{H}$ is to use Definition \ref{def:geomech} and add double-geometric noise with scale $2/\epsilon$ to each cell of $H$. However, because the maximum group size is not known, the length of $H$ is private. Thus we must determine a maximum non-private size $K$ and make the following modifications to $H$. If all groups in $H$ have less than $K$ people, then $H$ is extended with 0's until its length is $K+1$. If some groups of $H$ have size more than $K$, we change the sizes of those groups to $K$. Call the resulting histogram $H^\prime$. We then add double-geometric noise with scale $2/\epsilon$ to each cell of $H^\prime$.
This strategy satisfies $\epsilon$-differential privacy because the sensitivity of $H^\prime$ is 2:
\begin{lemma}
The global sensitivity of $H^\prime$ is 2.
\end{lemma}
\begin{proof}
For any group size $i<K$, adding a person to a group of size $i$ decreases $H^\prime[i]$ by one and increases $H^\prime[i+1]$ by one. When $i\geq K$, there is no change to $H^\prime$. Similarly, removing a person from a non-empty group of size $i\leq K$ increases $H^\prime[i-1]$ by one and decreases $H^\prime[i]$ by one (for a total change of 2). If $i>K$,  there is no change to $H^\prime$. 
\end{proof}

Let $\widetilde{H}$ be the noisy version of $H^\prime$. The numbers in $\widetilde{H}$ can be negative. Thus we post-process $\widetilde{H}$ to obtain an estimate $\widehat{H}$ by solving the following problem (using a quadratic program solver):
\begin{align*}
\widehat{H} &=\arg\min_{\widehat{H}} ||\widetilde{H} - \widehat{H}||_2^2\\
&\text{ s.t. }\widehat{H}[i]\geq 0 \text{ for all i}\quad \text{ and }\sum_i \widehat{H}[i] = G
\end{align*}
To get integers, we set $r = G-\sum_i \lfloor \widehat{H}[i] \rfloor$, round the cells with the $r$ largest fractional parts up, and round the rest down.


\changes{This approach has several weakness. First, there are many indexes $i$ where $H[i]=0$ and, after noise addition, many of them (a $\frac{2e^{-\epsilon}}{e^{-\epsilon}+1}$ fraction of indices) will have non-zero entries. Then roughly half of them end up with positive counts due to nonnegativity constraints. As a result, in our experiments, this method had several orders of magnitude worse error than the algorithms we describe next. Second, earthmover's distance is equivalent to measuring error between the cumulative sums of the estimate $\widehat{H}$ and the original data $H$: $\sum_{i} | \sum_{j\leq i} \widehat{H}[j] - \sum_{j\leq i} H[j]|$. Since noise is added to every cell independently, the $i^\text{th}$ component of the error depends on the sum of $i$ random variables (whose total variance is $O(i)$). Hence, if $H$ has $n$ components, one would expect the error to be $O(n^2)$.}

\subsection{Unattributed Histograms $H_g$}\label{subsec:unatt}

The next approach is to use algorithms for unattributed histograms \cite{Hay2010Boosting}. We can convert $H$ into the representation $H_g$, where $H_g[i]$ is the size of the $i^\text{th}$ smallest group. The length of $H_g$ is $G$, so potentially this can be a very large histogram. One of the properties of $H_g$ is that it is non-decreasing. Hence we can achieve $\epsilon$-differential privacy by adding independent double-geometric noise with scale $1/\epsilon$ to each element of $H_g$ to obtain $\widetilde{H}_g$, because the sensitivity of unattributed histograms is 1 \cite{Hay2010Boosting}. However, $\widetilde{H}_g$ is no longer non-decreasing (or even nonnegative), hence, following \cite{Hay2010Boosting,Lin2013Information}, we post-process it by solving the following optimization problem with either $p=1$ or $p=2$:
\begin{align*}
\widehat{H}_g &=\arg\min_{\widehat{H}_g} ||\widetilde{H}_g-\widehat{H}_g||_p^p\\
&\text{ s.t. } 0\leq \widehat{H}_g[i]\leq \widehat{H}_g[i+1] \text{ for }i=0,\dots, G-1
\end{align*}
Then we round each entry of $\widehat{H}_g$ to the nearest integer. Since the result was non-decreasing before rounding, it will remain non-decreasing after rounding.  Then from the resulting $\widehat{H}_g$ estimate, we convert it back to $\widehat{H}$ by counting, for each $i$, how many estimated groups have size $i$.

This optimization problem is known as \emph{isotonic regression}. When $p=2$ it can be solved in linear time using the min-max algorithm \cite{barlow1972isotonic}, pool-adjacent violators (PAV) \cite{barlowstatistical,robertson1968estimating}, or a commercial optimizer such as Gurobi \cite{gurobi}. When $p=1$ it can be converted to a linear program but it runs much slower. In our experiments,
we used $p=2$ because $H_g$ can have length in the hundreds of millions. For these sizes, the quadratic program is much faster to solve using PAV.


One observation we had is that this method is very good at estimating large group sizes, and most of its error comes from estimation errors of small group sizes (see Figure \ref{fig:abs_err}).

\begin{comment}
Formally, given a noisy sorted group array $\widetilde{H_g}$ and a set of constraints $\gamma$ where the total number of groups is enforced through the length of the array, we aim to find a sorted array  $\widehat{H_g}$ that is closest to $\widetilde{H_g}$ in terms of $L_2$ distance. The solution to this constrained $L_2$ minimization is equivalent to solving the isotonic regression with same weight on all data point.

\begin{definition}[$L_2$ Isotonic Regression on $H_g$] 
Given $\widetilde{H_g}$, the objective is to  find the solution $\widehat{H_g}$ which minimizes  ${|| \widetilde{H_g} - \widehat{H_g} ||_2}^2$ and subject to $\widehat{H_g}[i] \leq \widehat{H_g}[i+1]$, $0 \leq i <n$ and min($\widehat{H_g}$) = min($H_g$) $\geq 0$.
\end{definition}

The solution to this optimization is unique and can be solved by several algorithms such as the min-max algorithm \cite{barlow1972isotonic}. Because the sorted group size $H_g$ usually grows very large, we turn to use a faster algorithm pool-adjacent violator algorithm (PAV) \cite{barlowstatistical,robertson1968estimating} which specially designed for isotonic regression problem. Similarly, we round the estimated size to get $\widehat{H_g}$.
\end{comment}



\subsection{Cumulative Histograms $H_c$}\label{subsec:cumu}
Another strategy is to use $H_c$, the cumulative sum of $H$. Since the cumulative sum is non-decreasing, we again can add noise and use isotonic regression. As in the naive case, we must determine a public upper bound\footnote{\changes{Recall $K$ is an upper bound on the maximum number of people in a group. This method is not very sensitive to K -- in the experiments we used $K=100,000$ on datasets where the largest group had around 10,000 people -- an order of magnitude difference and still the estimated size of the largest group ended up being around 10,000. Thus if we have no prior knowledge, we can estimate $K$ as follows. Set aside a small privacy budget, since $K$ does not need much accuracy (e.g., $\epsilon=10^{-4}$). Let $X$ be the number of people in the largest group. Estimate $K$ as $X+$Laplace$(1/\epsilon) + 5\frac{\sqrt{2}}{\epsilon}$ -- i.e. add 5 standard deviations to a noisy estimate of $X$ so that $P(K \geq X)> 0.9995$ }} $K$ on $H$ and modify $H$ appropriately (as in Section \ref{subsec:naive}) before computing the cumulative sum. We saw in Section \ref{subsec:error} that error in estimation of \changes{count-of-counts\xspace histograms} is measured as the $L_1$ difference between cumulative size histograms. Thus it makes sense to privatize these histograms directly.
\begin{lemma}\label{lem:hcsens}
The global sensitivity of $H_c$ is 1.
\end{lemma}
\begin{proof}
Adding one person to a group of size $i$ means there is one less group of size $i$ and one more group of size $i+1$. Thus $H_c[i]$ decreases by one but $H_c[i+1]$ remains the same (i.e., number of groups $\leq i+1$ does not change). None of the other entries of $H_c$ change. Similarly, removing a person from a group of size $i$ means $H_c[i-1]$ increases by one and $H_c[i]$ does not change (and neither do any other entries). Thus the overall change is $1$.
\end{proof}
Thus we can satisfy differential privacy by adding independent double-geometric noise with scale $1/\epsilon$ to each cell of $H_c$ to get $\widetilde{H}_c$. We then postprocess $\widetilde{H}_c$ by solving the following optimization problem (using $p=1$ or $p=2$):
\begin{align*}
\widehat{H}_c &=\arg\min_{\widehat{H}_c} ||\widehat{H}_c - \widetilde{H}_c||_p^p\\
&\text{ s.t. } 0 \leq \widehat{H}_c[i] \leq \widehat{H}_c[i+1] \text{ for }i=0,\dots, K\\
& \text{ and }\widehat{H}[K]=G
\end{align*}
These problems can again be solved using PAV (in the case of $p=2$) or with a commercial optimizer (in the cases of $p=1$ or $p=2$).
In our experiments, we found that the $L_1$ version of the problem (with $p=1$) performs better than the $L_2$ version (with $p=2$). This is consistent with prior observations on unattributed histograms \cite{Lin2013Information}. A Bayesian post-processing is known to further reduce error, but we did not use it because it scales quadratically with the size of the histogram (the sizes of our histograms make this prohibitively expensive) \cite{Lin2013Information}. 
Finally, we then round $\widehat{H_c}$ to nearest integer and convert it back into a count-of-counts\xspace histogram $\widehat{H}$. We found that the $L_1$ version of the problem mostly returns integers, so rounding is minimal.

Unlike the $H_g$ method, we found that this method is accurate for small group sizes but less accurate for large groups.

\begin{figure}[t]
	\centering
	\includegraphics[width=7.5cm]{HgGood.png}
\vspace{-9pt}
	\caption{Error visualization. $x$-axis: cumulative sum of group size counts. $y$-axis: estimation error at these sizes. Top: $H_g$ method (Section \ref{subsec:unatt}). Bottom $H_c$ method (\ref{subsec:cumu}). Errors from the $H_g$ method are concentrated around small group sizes. Errors from $H_c$ are more dense throughout the rest of the sizes.}
    \label{fig:abs_err}
\end{figure}


\begin{comment}
\subsection{Dynamic Switching}\label{sec:switch}
Overall, the $H_c$ method provides more accurate estimates for smaller group sizes and much of its error comes from larger error for large group sizes. On the other hand $H_g$ method does the opposite -- it is more accurate for the larger groups and less accurate for the smaller groups. Since it is difficult to tell in advance which is more accurate, we employ a dynamic switching technique that is a variant of the parameter tuning algorithm proposed by Chaudhuri et al. \cite{Chaudhuri2011ERMJMLR}.

Specifically, given a privacy loss budget of $\epsilon$, we split it into $\epsilon_1$ and $\epsilon_2$ (so $\epsilon_1+\epsilon_2=\epsilon$). We then perform the following steps.
\begin{enumerate}
\item Use $\epsilon_1$ privacy budget to compute $\widehat{H}_g$ and convert it into a group size histogram $\widetilde{H}_1$.
\item Use $\epsilon_1$ privacy budget to compute $\widehat{H}_c$ and convert it into a group size histogram $\widetilde{H}_2$.
\item Using Earthmover's distance $d$, set $A = d(\widehat{H}_1, H)$ and $B=d(\widehat{H}_2, H)$. Compute $Z= A-B + X$ where $X$ is drawn from a double-geometric distribution with scale $2/\epsilon_2$.
\item If $Z > 0$, return ($\widehat{H}_1$, ``method $H_g$'') otherwise return ($\widehat{H}_2$, ``method $H_c$'').
\end{enumerate}
The difference between this algorithm and the tuning algorithm of Chaudhuri et al. is that we use double-geometric noise in the last step rather than the exponential mechanism.\footnote{It is part of folklore that using the exponential mechanism to select between $k$ items is equivalent to adding independent random variables from the standard Gumbel distribution to the score of each item, and picking the item with the largest noisy score. Our variant uses double-geometric noise rather than Gumbel and thus has less variance.}
A naive analysis would suggest that this procedure satisfies $(2\epsilon_1 + \epsilon_2)$-differential privacy. However, essentially the same proof as in \cite{Chaudhuri2011ERMJMLR} shows that it actually satisfies $\epsilon_1+\epsilon_2$-differential privacy.

We choose $\epsilon_2$ as follows. When $\epsilon_2=0.05$, then the variance of the Laplace distribution with scale $2/\epsilon_2$ has standard deviation $2\sqrt{2}/0.05\approx 56.6$ so most of the time we choose the best method $\widehat{H}_1$ or $\widehat{H}_2$ if their true errors are more then 3 standard deviations apart. If they are closer together, then both methods are comparable and either one is a good choice. However, when $\epsilon$ is small, we only want to use a small part of it to choose between $\widehat{H}_1$ and $\widehat{H}_2$. Thus we set $\epsilon_2 = \min(0.05, 0.05\epsilon)$ and $\epsilon_1 = \epsilon-\epsilon_2$.

\end{comment}

\subsection{Error Measure}\label{subsec:error}

The error measure is an important aspect of the problem as we would like to quantify the ``distance'' between $\tau.H$ and $\tau.\widehat{H}$. 
Ideally, we would like to measure this distance as the minimum number of people that must be added or removed from groups in $\tau.H$ to get $\tau.\widehat{H}$.

The standard error measures of Manhattan distance $||\tau.H - \tau.\widehat{H}||_1$ or sum-squared error $||\tau.H - \tau.\widehat{H}||_2^2$ do not capture this distance measure. To see why, suppose $H=[0, 100, 0,$ $ 0, 0, 0]$ meaning that all 100 groups have size $1$. Consider two estimates $\widehat{H}_1=[0, 0, 100, 0, 0, 0, 0]$ (where all groups have size 2) and $\widehat{H}_2=[0, 0, 0, 0, 0, 100]$ (where all groups have size 5). We see that $||H - \widehat{H}_1||_1=||H-\widehat{H}_2||_1=200$ and $||H - \widehat{H}_1||_2^2=||H-\widehat{H}_2||_2^2=20,000$. However,  we should consider $\widehat{H}_1$ to be closer to $H$ than $\widehat{H}_2$ is. The reason is that if we add one extra person into each group in $H_1$, we would obtain $\widehat{H}_1$. On the other hand, to obtain $\widehat{H}_2$, we would need to add 4 people to each group.

Thus, the appropriate way to measure distance between $H$ and $\widehat{H}$ is the Earthmover's distance (emd) as it precisely captures the number of people that must be added to or removed from groups in $H$ in order to obtain $\widehat{H}$. Normally, computing emd is not linear in the size of an array \cite{emd}. However, it can be computed in linear time for our problem by using the cumulative histograms.

\begin{lemma}[\cite{tcloseness}]\label{lem:error}
The earthmover's distance between $H$ and $\widehat{H}$ can be computed as $||H_c - \widehat{H}_c||_1$, where $H_c$ (resp., $\widehat{H}_c$) is the cumulative histogram of $H$ (resp., $\widehat{H}$). It is the same as the $L_1$ norm in the $H_g$ representation when the number of groups is fixed.
\end{lemma}

Our algorithms optimize error according to this metric. 


\subsection{Privacy Primitives}
We now describe the privacy primitives which serve as building blocks of our algorithm. One important concept is the sensitivity of a query.
\begin{definition}
Given a query $q$ (which outputs a vector), the global sensitivity of $q$, denoted by $\Delta(q)$ is defined as:
$$\Delta(q)=\max_{D_1,D_2}||q(D_1)-q(D_2)||_1,$$
where the maximum is taken over all databases $D_1, D_2$ that contain the public Hierarchy and Groups tables, and differ by the presence or absence of one record in the Entities table. 
\end{definition}

The sensitivity is used to calibrate the scale of noise needed to achieve differential privacy. We can use the Geometric Mechanism \cite{Ghosh2009Universal} instead of the Laplace Mechanism \cite{DMNS06Calibrating}, because we want our final counts to be integers. The Geometric mechanism is also preferable to the Laplace mechanism as it has lower variance, and is not susceptible to side-channel attacks when implemented in floating point arithmetic \cite{mironov12}.

\begin{definition}[Geometric Mechanism \cite{Ghosh2009Universal}]\label{def:geomech}

Given a \\ 
database $D$, a query $q$ that outputs a vector, and a privacy loss budget $\epsilon$, the geometric mechanism adds independent noise to each component of $q(D)$ using the following distribution: $P(X=k) = \frac{1-e^{-\epsilon}}{1+e^{-\epsilon}} e^{-\epsilon |k|/\Delta(q)}$ (for $k=0, \pm 1, \pm 2,$ etc.). This distribution is known as the double-geometric with scale $\Delta(q)/\epsilon$.
\end{definition}

\begin{lemma}[\cite{DMNS06Calibrating,Ghosh2009Universal}] The Geometric Mechanism satisfies $\epsilon$-differential privacy.
\end{lemma}

\begin{comment}

\begin{table}[htb] 
\center
\caption{Notations}
\label{tab:notation}
\begin{tabular}{ | c | c |    }
\hline
    Symbol     & Description    \\
 \hline
 $\Gamma$ & A k-ary tree\\
\hline
 $\tau$ & A tree node\\
 \hline
$H$ & Group size histogram \\
\hline
$H_c$ &  Cumulative sum of H    \\
\hline
$H_g$ & Sorted group size array\\
\hline
\\[-1em]
$\tilde{H}$, $\tilde{H}_c$, $\tilde{H}_g$ &  the noisy version of $H$, $H_c$, $H_g$\\
\hline
\\[-1em]
$\hat{H}$, $\hat{H}_c$, $\hat{H}_g$ &  the inference result for $H$, $H_c$, $H_g$\\
\hline
 $\gamma$ & min size, non-negativity, total constraints \\
\hline
$\beta$ & consistency constraint\\
\hline
\end{tabular}
\end{table}

A  k-ary tree histogram database $\Gamma$ is defined as a collection of nodes with $L$ levels where each node $\tau$ contains a group size histogram $H$.  A group can be associated with a household, a geographical region, etc. The size of a group $g$ is the number of instance in $g$. For a group size histogram $H$, $H[i]$ is the number of groups of size $i$. We denote $H_c$ as the cumulative sum of $H$ where $H_c[i]$ is the number of groups of size $\leq i$. The sorted group size array $H_g$ is a transformation of $H$ where all the groups are sorted by increasing size. $H_g[0]$ is the size of the smallest group, $H_g[1]$ is the size of the 2nd smallest group, etc. The length of the $H_g$ array equals the sum of the entries in $H$ (i.e., the total number of groups in $H$). For any of these representations $x$ (i.e., $H$, $H_c$, $H_g$), $\tilde{x}$ denotes $x$ with noise added independently to every entries.

We enforce constraints $\gamma$ on the data in every node $\tau$ where $\gamma$ includes the minimum size of the group, the non-negativity for each size, and the total number of groups in $t$. The consistency constraint $\beta$ in the application of a k-ary tree structure data is defined as the histogram at parent node is equal to the sum of histograms in all its children node. All notations are summarized in  \Cref{tab:notation}.



\subsection{Problem Definition}

Let $Q_i$ be the histogram of node $i$. The histogram in a parent node  $Q_p = \sum_j Q_p^j$ where $j$ is its child node.  For every node $p$ (except for a leaf node) in a k-ary tree, the consistency constraint $\beta$ requires $\hat{Q_p} = \sum_j \hat{Q_p^j}$. 

\label{sec:problem}
\begin{problem}[A k-ary Tree Consistent Estimation]
Suppose we have a k-ary tree, $\Gamma$ with $n$ nodes. Given a noisy k-ary tree $\tilde{\Gamma}$ where data in each node is noisy, the goal of post-processing is to find a new tree $\hat{\Gamma}$ with $n$ nodes that is closest to $\tilde{\Gamma}$ and, at the same time, satisfies the consistency constraints $\gamma$ and $\beta$. 
\end{problem}
\end{comment}

