# DP-SIPS: A simpler, more scalable mechanism for differentially private
partition selection

Marika Swanberg  Boston UniversityDepartment of Computer Science
<EMAIL_ADDRESS>,  Damien Desfontaines  Tumult Labs<EMAIL_ADDRESS>and
Samuel Haney  Tumult Labs<EMAIL_ADDRESS>
###### Abstract.

Partition selection, or set union, is an important primitive in differentially
private mechanism design: in a database where each user contributes a list of
items, the goal is to publish as many of these items as possible under
differential privacy.

In this work, we present a novel mechanism for differentially private
partition selection. This mechanism, which we call DP-SIPS, is very simple: it
consists of iterating the naive algorithm over the data set multiple times,
removing the released partitions from the data set while increasing the
privacy budget at each step. This approach preserves the scalability benefits
of the naive mechanism, yet its utility compares favorably to more complex
approaches developed in prior work.

Differential privacy, partition selection, scalable algorithms

††journalyear: YYYY††journalvolume: YYYY††journalnumber: X††doi:
XXXXXXX.XXXXXXX

##  1\. Introduction

Group-level aggregation is a fundamental building block for many data analysis
tasks. For example, doctors may want to understand mortality rates for
patients with different underlying conditions, and search engine developers
may want to study the frequency of different search queries that users are
making.

With increased awareness of the risks of data disclosure, increasingly data
analysts are using differentially private mechanisms to answer queries on
sensitive data sets. In both of the given examples, the aggregate values
(e.g., the mortality rates, or the frequencies of search queries) are
sensitive, but also the groups themselves are sensitive. Revealing that _some_
user in the data set has a particular disease or made a particular search
query could _itself_ violate privacy—even if it is never revealed _who_ that
person is. Additionally, the set of groups may be _a priori_ unknown or
impractical to enumerate, like in vocabulary extraction
(gopi2020differentially, ) or the private release of search queries
(korolova2009releasing, ).

In order to privately do a group-level aggregation in such settings, one must
first privately compute the set of groups represented in the data set—ideally,
releasing as many as possible while still maintaining privacy. More formally,
the mechanism $M$ gets a data set $x=(W_{1},\ldots,W_{N})$ where each user $i$
has a list of items $W_{i}$ and $M$ differentially privately releases a subset
$S\subseteq\cup_{i}W_{i}$ of the partitions that is as large as possible.

We study this problem, called _differentially private partition selection_
(also called key selection or set union), with a focus on approaches that can
be incorporated into general-purpose differentially private tooling. When
designing such frameworks, scalability is of paramount importance
(wilson2020differentially, ; amin2022plume, ; tumultanalyticswhitepaper, ):
underlying mechanisms must be able to run even when the input or the output is
too large to fit on a single machine. In such cases, the computation itself
must also be parallelizable, so it can run across multiple machines and avoid
unacceptably long running times.

Differentially private partition selection mechanisms have been proposed as
early as 2009 (korolova2009releasing, ), but recent work has shed a new light
on this problem, and proposed alternative approaches that bring significant
utility gains (gopi2020differentially, ; carvalho2022incorporating, ). To
obtain these utility improvements, these newer mechanisms use a _greedy_
approach: each user considers what items have been contributed by previous
users so far, and “chooses” which items to contribute according to a _policy_
, chosen carefully to maintain a sensitivity bound. Unfortunately, these
mechanisms do not scale: each user chooses their contribution based on the
contribution of all previous users, so the data has to be processed one user
after another, and the overall algorithm cannot be parallelized. This
intuition can be verified experimentally: we show that such algorithms
eventually time out or run out of memory when the input is very large, even on
clusters of multiple machines. Because these greedy algorithms do not scale,
differential privacy tools that need to handle large datasets cannot use these
smarter approaches, and instead must rely on the naive algorithm and its
underwhelming utility (privacyonbeam, ; pipelinedp, ; tumultanalyticssoftware,
).

This raises a natural question: can we achieve the utility benefits of policy-
based approaches, while preserving the scalability of more naive approaches?
In this work, we introduce a new approach that combines both benefits: DP-
SIPS, short for _scalable, iterative partition selection_.

DP-SIPS relies on a simple idea: rather than having to process the data of
each user sequentially, it runs the naive, massively-parallelizable algorithm
multiple times, splitting the privacy budget between each step and removing
items that were previously discovered. It uses a minuscule amount of privacy
budget to publish and then remove the “heavy hitters”–those which have an
enormous amount of weight in the histogram–and allocates the majority of the
budget to discovering items that remain after heavy hitters are removed from
the dataset. On skewed datasets, the removal of heavy hitters allows users to
allocate proportionally more weight in subsequent histograms to less-frequent
items while the increased privacy budget on later iterations also _lowers_ the
threshold for releasing an item.

As we show experimentally on multiple real-world and synthetic datasets, this
mechanism has a similar utility to greedy approaches, but scales horizontally:
increasing the number of cluster nodes significantly reduces runtime. This
makes it a suitable choice for implementation in general-purpose DP
infrastructure with high scalability requirements.

The rest of this paper is organized as follows:

  * •

In Section 2, we formally define the problem and the building blocks we use
for DP-SIPS and its privacy accounting.

  * •

In Section 3, we detail existing approaches to differentially private
partition selection.

  * •

In Section 4, we introduce our algorithm and the proof of its privacy
guarantees.

  * •

In Section 5, we report on the experimental evaluation of DP-SIPS.

  * •

Finally, in Section 6, we discuss our results, report on some unsuccessful
approaches that we tried, and outline directions for future work.

##  2\. Preliminaries

A data set $x=(W_{1},\ldots,W_{N})$ contains a set of user lists
$W_{i}\in\mathcal{U}^{*}$. We refer to elements in $\mathcal{U}$ as _items_ or
_partitions_ , and define partition selection (also called key selection or
set union) as follows:

###### Definition 1 (Private Partition Selection Problem).

Given a (possibly unbounded) universe $\mathcal{U}$ of items and a data set
$x=(W_{1},\ldots,W_{N})$ of user lists $W_{i}\in\mathcal{U}^{*}$, an algorithm
$\mathcal{M}$ solves the _private partition selection problem_ if it is
differentially private, and $\mathcal{M}$ outputs a set
$S\subseteq\cup_{i}W_{i}$.

We begin by presenting the standard notion of differential privacy. Two data
sets $x,x^{\prime}$ are _neighbors_ if they differ on one user’s list:
$x=x^{\prime}\cup W_{i^{*}}$. Informally, differential privacy requires that
an algorithm’s output is distributed similarly on every pair of neighboring
data sets.

###### Definition 2 (Differential Privacy (dwork2006our, ;
dwork2006calibrating, )).

A randomized algorithm $\mathcal{M}:\mathcal{U}^{*}\to\mathcal{Y}$ is
$(\varepsilon,\delta)$-differentially private if for every pair of neighboring
datasets $x,x^{\prime}\in\mathcal{U}^{*}$ and for all subsets
$Y\subseteq\mathcal{Y}$,

$\Pr[\mathcal{M}(x)\in Y]\leq
e^{\varepsilon}\cdot\Pr[\mathcal{M}(x^{\prime})\in Y]+\delta.$

A common variant of differential privacy, called zCDP, is useful for analyzing
algorithms that sample noise from a Gaussian distribution (as ours will). The
definition of zCDP uses the Rényi Divergence:

###### Definition 3 (Rényi Divergence).

Fix two probability distributions $P$ and $Q$ over a discrete domain $S$.
Given a positive $\alpha\neq 1$, Rényi divergence of order $\alpha$ of
distributions $P$ and $Q$ is

$D_{\alpha}(P||Q)=\frac{1}{1-\alpha}\log\left(\sum_{s\in
S}P(s)^{\alpha}Q(s)^{1-\alpha}\right).$

###### Definition 4 ($\rho$-zCDP (bun2016concentrated, )).

A randomized mechanism $\mathcal{M}:\mathcal{X}^{*}\to\mathcal{Y}$ satisfies
$\rho$-zCDP if, for all $x,x^{\prime}\in\mathcal{X}^{*}$ differing on a single
entry,

(1) $\displaystyle
D_{\alpha}(M(x)||M(x^{\prime}))\leq\rho\cdot\alpha\quad\forall\alpha\in(1,\infty).$

This definition can also be relaxed to _approximate_ zCDP.

###### Definition 5 (Approximate zCDP (bun2016concentrated, )).

A randomized mechanism $M:\mathcal{X}^{n}\to\mathcal{Y}$ is
$\delta$-approximately $\rho$-zCDP if, for all
$x,x^{\prime}\in\mathcal{X}^{n}$ differing on a single entry, there exist
events $E=E(M(x))$ and $E^{\prime}=E^{\prime}(M(x^{\prime}))$ such that, for
all $\alpha\in(1,\infty)$,

(2) $\displaystyle
D_{\alpha}(M(x)|_{E}||M(x^{\prime})|_{E^{\prime}})\leq\rho\cdot\alpha\quad\text{and}$
(3) $\displaystyle
D_{\alpha}(M(x^{\prime})|_{E^{\prime}}||M(x)|_{E})\leq\rho\cdot\alpha,$

and $\Pr[E]\geq 1-\delta$ and $\Pr[E^{\prime}]\geq 1-\delta$.

Approximate zCDP satisfies composition and post-processing properties.

###### Lemma 6 ((bun2016concentrated, ), Lemma 8.2).

Let $M:\mathcal{X}^{n}\to\mathcal{Y}$ and
$M^{\prime}:\mathcal{X}^{n}\times\mathcal{Y}\to\mathcal{Z}$. Suppose $M$
satisfies $\delta$-approximate $\rho$-zCDP and for all $y\in\mathcal{Y}$,
$M^{\prime}(\cdot,y):\mathcal{X}^{n}\to\mathcal{Z}$ satisfies
$\delta^{\prime}$-approxmiate $\rho^{\prime}$-zCDP. Define
$M^{\prime\prime}:\mathcal{X}^{n}\to\mathcal{Z}$ by
$M^{\prime\prime}(x)=M^{\prime}(x,M(x))$. Then, $M^{\prime\prime}$ satisfies
$(\delta+\delta^{\prime}-\delta\cdot\delta^{\prime})$-approximate
$(\rho+\rho^{\prime})$-zCDP.

Next we show a conversion from approximate zCDP to approximate DP. We rewrite
the conversion lemma from (bun2016concentrated, ) to be slightly more general
(i.e. it uses an existing conversion from zCDP to approximate DP), and then
apply the tight zCDP-to-approximate DP conversion given in
(canonne2020discrete, ). Overall, this gives a tighter approximate zCDP-to-
approximate DP conversion than what is stated in Lemma 8.8 of
(bun2016concentrated, ).

###### Lemma 7 (Generalized version of Lemma 8.8 from (bun2016concentrated,
)).

Suppose we can show that every mechanism that satisfies $\rho$-zCDP must
satisfy $\epsilon^{*}(\rho),\delta^{*}(\rho)$-approximate DP. That is,
($\epsilon^{*},\delta^{*}$) is a function converting a (pure) zCDP guarantee
to an approximate DP guarantee. Suppose
$\mathcal{M}:\mathcal{X}^{N}\to\mathcal{Y}$ satisfies $\delta$-approximate
$\rho$-zCDP. Then, $\mathcal{M}$ satisfies
$(\epsilon^{*}(\rho),\delta+(1-\delta)\delta^{*}(\rho))$-DP.

###### Lemma 8 ((canonne2020discrete, ), Proposition 7).

Suppose $\mathcal{M}:\mathcal{X}^{N}\to\mathcal{Y}$ satisfies $\rho$-zCDP.
Then $\mathcal{M}$ satisfies $(\epsilon,\delta)$-approximate DP for any
$\epsilon>0$ and

(4)
$\delta=\inf_{\alpha\in(1,\infty)}\frac{\exp((\alpha-1)(\alpha\cdot\rho-\epsilon))}{\alpha-1}\left(1-\frac{1}{\alpha}\right)^{\alpha}.$

Combining Lemma 7 and Lemma 8 gives us the following.

###### Corollary 9.

Suppose $\mathcal{M}:\mathcal{X}^{N}\to\mathcal{Y}$ satisfies
$\delta$-approximate $\rho$-zCDP. Then $\mathcal{M}$ satisfies
$(\epsilon,\delta+(1-\delta)\delta^{\prime})$-approximate DP for any
$\epsilon>0$ and

(5)
$\delta^{\prime}=\inf_{\alpha\in(1,\infty)}\frac{\exp((\alpha-1)(\alpha\cdot\rho-\epsilon))}{\alpha-1}\left(1-\frac{1}{\alpha}\right)^{\alpha}.$

A common primitive in building private algorithms, the Gaussian Mechanism,
satisfies $\rho$-zCDP.

###### Definition 10 (Gaussian Distribution).

The Gaussian distribution with parameter $\sigma$ and mean 0, denoted
$\mathcal{N}(0,\sigma^{2})$ is defined for all $\ell\in\mathbb{R}$ and has
probability density

$h(\ell)=\frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{\ell^{2}}{2\sigma^{2}}}.$

###### Definition 11 ($\ell_{2}$-Sensitivity).

Let $f:\mathcal{U}^{n}\to\mathbb{R}^{d}$ be a function. Its
$\ell_{2}$-sensitivity is

$\Delta_{f}=\max_{\begin{subarray}{c}x,x^{\prime}\in\mathcal{U}\\\
x,x^{\prime}neighbors\end{subarray}}\|f(x)-f(x^{\prime})\|_{2}.$

###### Definition 12 (Gaussian Mechanism (bun2016concentrated, )).

Let $f:\mathcal{U}^{n}\to\mathbb{R}^{d}$ be a function with
$\ell_{2}$-sensitivity $\Delta_{f}$. Then the Gaussian mechanism is the
algorithm

$\mathcal{M}_{f}(x)=f(x)+(Z_{1},\ldots,Z_{d}),$

where $Z_{i}\sim\mathcal{N}\left(0,\frac{\Delta_{f}^{2}}{2\rho}\right)$.
Algorithm $\mathcal{M}_{f}$ satisfies $\rho$-zCDP.

##  3\. Prior Approaches to Partition Selection

In this section we discuss three existing algorithms for differentially
private partition selection. We begin with the naive algorithm, called
Weighted Gaussian, in Section 3.1. In Section 3.2, we present Policy Gaussian
(gopi2020differentially, ) and Greedy updates Without sampling
(carvalho2022incorporating, ).

The algorithms all have three main steps: first, they compute a _weighted
histogram_ , which is simply a mapping from an item $u\in\mathcal{U}$ to a
weight $H[u]\in\mathbb{R}$; second, they add calibrated noise to each item in
the histogram; and lastly, they release items that are above some
appropriately-chosen threshold $T$. The primary difference between the
algorithms is in how they compute the weighted histogram. This high-level
algorithm for private partition selection is described in Algorithm 1, which
can be composed with different weighted histogram algorithms.

Algorithm 1 High-Level Partition Selection Algorithm

###  3.1. Baseline: Weighted Gaussian

The Weighted Gaussian algorithm (Algorithm 2) is one of the simplest and first
algorithms for private partition selection (korolova2009releasing, ). To build
a weighted histogram, the algorithm first pre-processes the data set by
removing any duplicates within a user’s set and truncating each user’s set to
have at most $\Delta_{0}$ items. Then, the users compute a histogram with
bounded $\ell_{2}$-sensitivity as follows: each user $i$ updates the weight
$H[u]$ for each of the items $u$ in their set $W_{i}$ with the following rule:

$H[u]\leftarrow H[u]+\frac{1}{\sqrt{|W_{i}|}}.$

The resulting weighted histogram has an $\ell_{2}$-sensitivity of 1, so it can
be composed with the high-level algorithm using calibrated Gaussian noise and
an appropriate threshold to ensure that the overall algorithm satisfies
$\delta$-approximate $\rho$-zCDP. In the statement of Algorithm 2,
$\Phi(\cdot)$ is used to denote the cumulative density function of the
standard Gaussian distribution and $\Phi^{-1}(\cdot)$ is its inverse.

Algorithm 2 Weighted Gaussian

###### Theorem 1.

Fix any $\rho>0$, any $\delta\in(0,1)$, and any $\Delta_{0}\in\mathbb{N}$. The
Weighted Gaussian algorithm (Algorithm 2) satisfies $\delta$-approximate
$\rho$-zCDP.

See Appendix A for the proof of Theorem 1.

The Weighted Gaussian algorithm benefits from being highly scalable. In
particular, it lends itself well to parallel computation across several
computers within a cluster, because the weights on each histogram item can be
computed in parallel as well as the noise addition and thresholding steps (see
Figure 8). Thus, Algorithm 2 is the standard approach for doing partition
selection on data sets that are too large to fit in a single machine’s memory.
Unfortunately, Weighted Gaussian suffers from poor accuracy compared to the
greedy approaches we discuss next.

###  3.2. Greedy Approaches

One problem with the Weighted Gaussian algorithm is users waste their
sensitivity budget on histogram items that are already well above the
threshold. Most real-world data has highly skewed item frequencies, but
Weighted Gaussian increments all items in a user’s set by the same amount.

The two greedy algorithms we discuss next solve this problem by iterating
through the users one-by-one and using an _update policy_ and the current
state of the histogram to decide how to allocate weight across the items in
their set. That is, each user’s update depends on previous users’ updates. For
example, in both algorithms, users do not contribute to items in $H$ that have
already reached $T^{*}$, the threshold $T$ plus some positive buffer; items
that have reached the buffered threshold are very likely to be returned after
noise is added and thus do not need more weight. These adaptive update rules
are carefully chosen so that the overall algorithm has bounded global
sensitivity. As we will discuss in later sections, the main downside of these
algorithms is their sequential nature. By design, they require iterating over
the entire data set, which may be prohibitively slow for industrial data sets.

####  3.2.1. Policy Gaussian (gopi2020differentially, ) (DPSU)

As with the Weighted Gaussian algorithm, the data set is pre-processed by
removing duplicate items from each user’s set and truncating the user sets to
some fixed maximum size $\Delta_{0}$. Then, the algorithm iterates
sequentially over the users and, for each item $u$ in user $i$’s set $W_{i}$,
the user increments $H[u]$ by a weight that is proportional to $T^{*}-H[u]$,
where $T^{*}$ is equal to $T$ plus some positive buffer. Essentially, _items
that are further from the buffered threshold get more weight added to them,
and those that have already reached the buffered threshold get none_. The
weight that a user adds to each item is normalized so a single user’s update
to $H$ has an $\ell_{2}$-norm of at most 1.We refer to this algorithm as DPSU.

Gopi et al. prove that the global $\ell_{2}$-sensitivity of the entire Policy
Gaussian algorithm is bounded by 1, so applying the high-level algorithm
(Algorithm 1) with appropriately-scaled Gaussian noise to each item and
thresholding are sufficient for satisfying differential privacy.

####  3.2.2. Greedy updates Without sampling (GW) (carvalho2022incorporating,
)

Carvalho et al. observe that, rather than removing duplicate items in a user’s
set, one can use this frequency information to decide where to allocate
sensitivity budget. GW iterates over the users, and each user computes
$u^{*}$: the most frequent item in their list such that $H[u^{*}]$ is below
the buffered threshold $T^{*}$. Then, the user increments $H[u^{*}]$ by
$\min(1,T^{*}-H[u^{*}],budget)$ where $budget$ is the user’s remaining
$\ell_{1}$ budget. This process is repeated until the user’s initial budget of
1 is consumed or the user has no items left that are below $T^{*}$111As the
algorithm is presented in (carvalho2022incorporating, ), it does not always
terminate. In particular, they do not consider the case when a user has budget
left but all of its items have reached $T^{*}$, but adding this condition
luckily does not affect the privacy proof.. Carvalho et al. prove that their
GW algorithm for building a weighted histogram has a global
$\ell_{1}$-sensitivity bound of 1, so running the high-level algorithm
(Algorithm 1) with Laplace noise and thresholding is sufficient for satisfying
differential privacy.

One notable feature of their algorithm is that it can use item frequency
information from a public data set to increase the accuracy of the frequency
estimates. We do not use this feature when comparing it to other algorithms,
since our goal is to create a general-purpose algorithm that could be
incorporated into a privacy framework without requiring a data analyst to
input a public data set.

##  4\. DP-SIPS

Figure 1. Depiction of Weighted Gaussian noisy histogram (left) compared to
intermediate SIPS noisy histograms (right three diagrams) on a skewed data
set. Solid blue bars represent partitions that will be returned, and yellow
outlines represent the weight of each item on the previous iteration. Although
the threshold for Weighted Gaussian is lower than each of the SIPS thresholds,
SIPS benefits from less-frequent items getting increased weight as returned
partitions are removed from user sets.

In this section we present DP-SIPS: Differentially Private Scalable, Iterative
Partition Selection, detailed in Algorithm 3. The basic structure is quite
simple: it runs Weighted Gaussian on the data set multiple times with
increasing privacy budget (and corresponding decreasing thresholds); on each
iteration, the partitions returned by Weighted Gaussian are removed from each
of the users’ sets.

Because Weighted Gaussian updates the histogram uniformly across a user’s
items, when returned partitions are removed from the user’s set, they can
allocate more weight to their remaining items (see Figure 1). The first
iteration returns the very popular items using only a tiny fraction of the
overall privacy budget, and subsequent iterations yield less frequent items.
The action of the algorithm is twofold: _in each iteration, the threshold is
lowered at the same time as users allocate more weight to each of the items
that remain in their sets_ (since the user sets get smaller when previously-
returned items are removed).

Furthermore, the user sets are re-truncated on each iteration after the
previously returned partitions are removed. For data sets that are both skewed
in the item frequencies and in the sizes of users’ sets, this allows
additional items to be included on each iteration.

Algorithm 3 DP-SIPS: Scalable, Iterative Partition Selection

###### Theorem 1.

For any $\rho>0$ and any $\delta\in(0,1]$, Algorithm 3 satisfies
$\delta$-approximate $\rho$-zCDP.

###### Proof.

By Theorem 1, each call to Weighted_Gauss satisfies $\delta_{i}$-approximate
$\rho_{i}$-zCDP. Applying composition and postprocessing (Lemma 6), Algorithm
3 satisfies $\left(\sum_{i=0}^{I-1}\delta_{i}\right)$-approximate
$\left(\sum_{i=0}^{I-1}\rho_{i}\right)$-zCDP.

If $r=1$, then clearly $\sum_{i=0}^{I-1}\rho/I=\rho$ and similarly
$\sum_{i=0}^{I-1}\delta/I=\delta$.

Now, for $r\neq 1$ we will solve for these summations using the closed-form
formula for geometric sums.

$\displaystyle\sum_{i=0}^{I-1}\rho_{i}$
$\displaystyle=\sum_{i=0}^{I-1}\rho\cdot
r^{I-i-1}\cdot\frac{1-r}{1-r^{I}}=\rho\cdot\frac{1-r}{1-r^{I}}\cdot\sum_{i=0}^{I-1}r^{I-i-1}$
$\displaystyle=\rho\cdot\frac{1-r}{1-r^{I}}\cdot\sum_{j=0}^{I-1}r^{j}=\rho\cdot\frac{1-r}{1-r^{I}}\cdot\frac{1-r^{I}}{1-r}=\rho.$

An identical calculation holds for $\sum_{i=0}^{I-1}\delta_{i}$. Thus,
Algorithm 3 satisfies $\delta$-approximate $\rho$-zCDP. ∎

A key advantage of DP-SIPS is that it is massively parallelizable: within each
iteration, the mechanism performs a simple groupby-count operation, which can
be distributed among multiple cores or a cluster of multiple machines. DP-SIPS
requires multiple iterations, so it takes longer to run than the naive
algorithm. However, only a few iterations are required to achieve good
accuracy (see Section 5.2.3), so its runtime is still reasonable even on large
datasets. A visual explanation of the scalability property of DP-SIPS in
comparison with other approaches can be found in Figure 8, and an experimental
performance and scalability comparison can be found in Section 5.3.

##### Relation to Algorithm 7 in (kamath2019privately, )

Private Product-Distribution Estimator in (kamath2019privately, ) privately
estimates a product of $d$ Bernoulli distributions, and has some structural
similarities to DP-SIPS: it also uses weighted Gaussian iteratively with
decreasing thresholds, removing items above the threshold in each round to
estimate the frequencies of each coordinate. Aside from the problem setting
(bounded vs. unbounded domain) and goal (distribution estimation vs. partition
selection), a major difference is that (kamath2019privately, ) partitions the
dataset horizontally into $\log(d/2)$ groups, and uses different subsets of
user records on each round. This partitioning step is key to their privacy and
accuracy arguments and leads to $\log(d/2)$ rounds rather than the constant
number of rounds that DP-SIPS uses.

###  4.1. Scalability Analysis in the MapReduce Framework

We analyze DP-SIPS in the MapReduce model of distributed computing to better
understand its scaling behavior on a cluster. The model proceeds in rounds
that are synchronized across worker nodes running in parallel. The data
records begin arbitrarily partitioned among the worker nodes. Each round
begins with a _map_ phase in which the worker node applies a map to each
record in its partition. Next, the records are _shuffled_ among the workers so
that records with the same key are located on the same worker node. Lastly,
each worker processes the records in the _reduce_ phase. The three phases
(map, shuffle, and reduce) constitute one round in the MapReduce framework.
Typically, the shuffle phase is the most time intensive since large amounts of
data must move among the worker nodes, so the overall number of rounds usually
reflects an algorithm’s efficiency on a real system.

DP-SIPS proceeds in the following steps (see Figure 2):

  1. (1)

Shuffle records by user_id

  2. (2)

Bound the number of contributions per user

  3. (3)

Shuffle records by item. Call this dataset $D$.

  4. (4)

Count the number of items, add noise, and store the set of items $S$ above the
threshold

  5. (5)

Do a join to remove items in $S$ from $D$. The records now consist of
$D\setminus S$. If this is the last iteration, this step can be skipped.

  6. (6)

Return to Step 1 for $I$ iterations in total.

All but the last iteration of DP-SIPS involves three MapReduce rounds: first
shuffling records by user to truncate; then shuffling by item to count, add
noise, and threshold; lastly, doing a join with the original data to remove
items that were released. In the last iteration, the join is not necessary.
So, in total the algorithm does $3(I-1)+2$ rounds in the MapReduce framework,
where $I$ is the number of iterations of DP-SIPS. In our experiments we set
$I=3$, so DP-SIPS performs 8 MapReduce rounds. Weighted Gaussian only uses 2
MapReduce rounds, since $I=1$. The greedy algorithms cannot be
parallelized–specifically the greedy methods for histogram computation–so each
record must be loaded into the one worker node to make its contribution to the
histogram. This requires $O(N)$ rounds of MapReduce, where $N$ is the number
of users. Indeed, our experimental results in Section 5.3 (Figures 10 and 11)
show that the runtime of DP-SIPS decreases with the number of cores while the
greedy algorithms stay constant.

Figure 2. MapReduce Diagram for DP-SIPS. Each iteration within SIPS involves
two rounds of MapReduce.

##  5\. Experimental Results

We use data sets of several different sizes from varied text domains to
validate the empirical performance of the DP-SIPS algorithm, and we find that
in general DP-SIPS has comparable accuracy to DPSU and GW while scaling to
large data sets that DPSU and GW timeout on. As with prior works, we focus on
the problem of vocabulary discovery under the constraint of user-level
privacy.

We begin by describing the data sets in Section 5.1 and then, in Section 5.2
we discuss how the empirical accuracy of DP-SIPS compares to that of existing
algorithms, and in Section 5.3 we discuss the results from our scalability
experiments.

###  5.1. Data sets

We use six publicly-available data sets and one synthetically-generated data
set to study the accuracy of DP-SIPS against existing partition selection
algorithms, and we use the largest two of the seven for scalability
experiments on Amazon Elastic Map Reduce clusters. Reddit (reddit_dataset, )
is a data set of text posts collected from r/AskReddit which appeared as a
benchmark in (gopi2020differentially, ; carvalho2022incorporating, ). We also
use four data sets that appeared in (carvalho2022incorporating, ): Twitter
(twitter_dataset, ), comprising customer support tweets to and from large
corporations; Finance (finance_dataset, ), financial headlines for stocks;
IMDb (imdb_dataset, ), a set of movie reviews scraped from IMDb; Wikipedia
(wiki_dataset, ), a set of Wikipedia abstracts (where we treat each abstract
as a separate user set).

For the scalability experiments, we use Amazon (amazon_dataset, ), a publicly-
available text data set from Kaggle of 4 million Amazon product reviews, and a
synthetically generated data set Synthetic_80M with 80 million users and 4.6
billion observations. The synthetic data set, Synthetic_80M, is generated as
follows. First, each user draws the number of items in their set from a Pareto
distribution of scale 10 and shape 1.16; this captures the common “80-20
principle” observed on many empirical datasets that states that 80% of
outcomes are due to 20% of causes. Then, each item in each user’s set is
generating by sampling from a zeta distribution of parameter $1.1$; this
captures Zip’s law, which states that in many types of data (including words
in natural languages), the rank-frequency distribution follows an inverse
relation.

Table 1 lists the number of users, number of observations, and the vocabulary
size of each data set. Using the same methods as (gopi2020differentially, ;
carvalho2022incorporating, ), we preprocess all non-synthetic data sets using
tokenization with nltk.word_tokenize, removing URLs and symbols, and lower
casing all words.

Data set | Users | Observations | Vocabulary Size
---|---|---|---
Reddit | 223,388 | 373,983 | 155,701
Twitter | 702,682 | 2,811,774 | 1,300,123
Finance | 1,400,465 | 1,400,465 | 267,256
IMDb | 49,999 | 49,999 | 194,532
Wikipedia | 245,103 | 245,103 | 631,866
Amazon | 4,000,000 | 4,000,000 | 4,250,427
Synthetic_80M | 80,000,000 | 4,643,596,660 | 741,129,124
Table 1. Number of users, number of observations, and true vocabulary size for
each data set we consider.

###  5.2. Accuracy results

Because the goal of partition selection is to privately output as many
partitions as possible, we measure accuracy as the _number of partitions
released_. To date, there are no analytical accuracy guarantees for any prior
partition selection algorithms (in the setting where users contribute multiple
items), so we must use experimental validation to understand the accuracy of
each algorithm.

We test the accuracy of the four algorithms: Weighted Gaussian (Algorithm 2),
SIPS (our Algorithm 3), Policy Gaussian (DPSU) from (gopi2020differentially,
), and Greedy updates Without sampling (GW) from (carvalho2022incorporating,
).

Table 2 shows the following accuracy trends for SIPS on the datasets described
in Section 5.1:

  * •

SIPS only performs slightly worse than both GW and DPSU on one dataset
(Finance); in general, SIPS’ performance is on par with DPSU’s.

  * •

The accuracy of SIPS is consistently approximately double that of Weighted
Gaussian.

Data set | Wt. Gauss | SIPS | DPSU | GW
---|---|---|---|---
Reddit | 6,160 | 11,392 | 11,186 | 11,984
Twitter | 12,632 | 23,649 | 23,576 | 27,184
Finance | 17,350 | 27,559 | 29,005 | 37,503
IMDb | 3,728 | 7,759 | 5,845 | 3,133
Wikipedia | 11,340 | 21,037 | 18,129 | 11,251
Amazon | 67,522 | 144,805 | 143,997 | 185,563
Synthetic_80M | 711,601 | 1,137,467 | - | -
Finance (deduplicated) | - | - | - | 37,563
IMDb (deduplicated) | - | - | - | 3,005
Wikipedia (deduplicated) | - | - | - | 9,802
Table 2. Number of partitions returned by Wt. Gauss (Algorithm 2), SIPS
(Algorithm 3), DPSU (Policy Gaussian from (gopi2020differentially, )), and GW
(from (carvalho2022incorporating, )) on six data sets. Additionally, we run GW
on three data sets where duplicates within user lists have been removed. Note
that the other three algorithms already deduplicate user sets. For SIPS we use
3 iterations and $r=1/3$. For Wt. Gaussian, SIPS, and DPSU, the privacy budget
is set to $\rho=0.1,\delta=10^{-5}$, and the user contributions are truncated
to $\Delta_{0}=100$. For GW, $\varepsilon=1.7$ and $\delta=8.1142\times
10^{-5}$, which implies $10^{-5}$-approximate 0.1-zCDP. On Synthetic_80M, DPSU
and GW ran out of memory before completing the computation.

In addition, Figures 3, 4, and 5 demonstrate that the relative accuracy of
each algorithm on a given data set remains consistent across different choices
of privacy parameters.

####  5.2.1. Accuracy Inconsistencies

While GW tends to perform well, its accuracy on the IMDb and Wikipedia data
sets is below even that of Weighted Gaussian (see Figure 5). To further
investigate this phenomenon, we modify these two data sets in addition to
Finance to remove repeated items within each user’s set (see the data sets
marked as “deduplicated” in Table 2). Note also that the accuracy of the other
three algorithms is unaffected by deduplication since they all perform this
preprocessing step to the data sets.

Without any item frequency information, the algorithm adds all of its weight
to a randomly-selected item, so one would expect the performance on all three
data sets to be worse than Weighted Gaussian. To the contrary, GW’s accuracy
on the deduplicated Finance data set is still significantly higher than the
others, and GW’s accuracy on the deduplicated IMDb and Wikipedia is again
worse than Weighted Gaussian. Figure 5 shows that GW’s poor relative accuracy
on IMDb worsens at higher levels of epsilon.

We believe that GW’s inconsistent accuracy depends on the ratio between
vocabulary size and number of observations: GW performs very well on datasets
where this ratio is small (Finance: 0.19), and very poorly on datasets with
large relative vocabulary (IMDb: 3.8). Other factors, such as the short length
of Finance headlines compared to IMDb reviews and Wikipedia abstracts, might
also play a role. We did not attempt to fully resolve this discrepancy in GW’s
accuracy; however, we present these results as a caution when selecting a
partition selection algorithm.

Figure 3. Accuracy on Reddit as a function of varying epsilon (and associated
$\rho$), for $\delta=10^{-5}$. Figure 4. Accuracy on Reddit as a function of
varying delta for $\varepsilon=1.7$. Figure 5. Accuracy on IMDb as a function
of varying epsilon for fixed $\delta=10^{-5}$.

####  5.2.2. Comparing zCDP to DP

We implement our algorithm to satisfy approximate zCDP to take advantage of
its simpler composition properties over standard DP. In
(gopi2020differentially, ), the Policy Gaussian algorithm satisfies
$(\varepsilon,\delta)$-DP, but the threshold and Gaussian noise can easily be
recalibrated to satisfy approximate zCDP instead. We do so in this work to
facilitate the comparison to our algorithm. Unfortunately, the GW algorithm
from (carvalho2022incorporating, ) has a bounded $\ell_{1}$-sensitivity and
uses Laplace noise, so it is not easily converted to satisfy approximate zCDP
without a significant loss in accuracy. Instead, we use Corollary 9 to choose
appropriate $(\varepsilon,\delta)$ parameters for the given $\rho$ and
$\delta_{CDP}$. The tables in Appendix B give the conversions derived from
Corollary 9 used for our experiments.

The conversion from approximate zCDP to approximate DP is not exactly tight
(meaning the given $(\varepsilon,\delta)$ may be higher than the true privacy
guarantee given by the approximate zCDP parameters). Because of the looseness
of the conversion, the GW algorithm’s accuracy may be slightly inflated when
compared to the other algorithms.

Doing the privacy budget accounting using approximate zCDP is not exactly
tight either; instead, one could use numerical methods to compute the total
privacy budget spent (meiser2018tight, ). At the time of writing, we could not
find a working implementation of such methods that could support the privacy
property of Weighted Gaussian. We do not believe that the results of the
experiments would be meaningfully different by doing the privacy analysis this
way. Furthermore, using approximate zCDP has the advantage of being easier to
integrate with existing differential privacy software
(tumultanalyticssoftware, ; tumultanalyticswhitepaper, ).

Doing the privacy budget accounting with zCDP is not exactly tight either:
numerical methods using

####  5.2.3. Selecting hyperparameters

Our algorithm has several hyperparameters (aside from the privacy parameters
$\rho$ and $\delta$) that need to be set by the data analyst: $\Delta_{0}$ the
maximum number of items per user, $r$ the ratio between the privacy budget for
iteration $i+1$ and $i$, and $I$ the number of iterations in the algorithm.
One option is to divide the privacy budget and try several hyperparameter
settings and select the setting with the highest accuracy; however, this
wastes a lot of privacy budget. We find that the accuracy of DP-SIPS is
largely invariant to reasonable settings of the hyperparameters, and we
provide general rules of thumb for selecting them.

Figures 6 and 7 display the accuracy of our algorithm on five data sets with
varying $r$ and $\Delta_{0}$, respectively. For the given setting of $\delta$
and $\rho$, the accuracy of DP-SIPS is largely unaffected by different choices
of $r$ and $\Delta_{0}$ across 5 data sets, and Table 3 shows that the
accuracy of DP-SIPS only slightly increases with the number of iterations $I$
on the Reddit data set (for the given settings of the other parameters).
Furthermore, Figure 6 (privacy ratio $r=1$) suggests that evenly splitting the
budget among the rounds yields lower accuracy than geometrically increasing
the budget at each round.

The results suggest that DP-SIPS has good accuracy for a large range of
hyperparameters (aside from the privacy parameters). We recommend using the
following settings: $r\in[0.2,0.4]$, $I\geq 3$, and $\Delta_{0}$ set to an
overestimate of the true maximum number of per-user contributions. If the true
maximum number of contributions per user is public, $\Delta_{0}$ should be set
to that value.

Iterations $I$ | Partitions Returned
---|---
1 | 5,464
2 | 10,041
3 | 11,126
4 | 11,182
5 | 11,541
6 | 11,061
8 | 11,585
10 | 11,637
Table 3. SIPS’s accuracy as a function of the number of iterations on Reddit
with parameters $\rho=0.1,\delta=10^{-5}$, $r=1/3$, and $\Delta_{0}=50$.
Figure 6. Number of partitions returned by SIPS versus the choice of the
privacy ratio parameter, $r$. All other hyperparameters are identical to those
listed in Table 2. For all 5 data sets, the accuracy is maximized for
$r\in[0.1,0.4]$. Figure 7. Number of partitions returned by SIPS for varying
per-user maximum contribution bounds, $\Delta_{0}$. All other hyperparameters
are identical to those listed in Table 2. On all 5 data sets, increasing
$\Delta_{0}$ past 100 has almost no affect on accuracy, with maximum accuracy
for $\Delta_{0}\in[5,150]$. We note that Finance is unusual as each financial
headline is quite short.

###  5.3. Scalability results

Figure 8. Parallelism diagram for Weighted Gaussian (left), SIPS (middle), and
greedy algorithms (right). Arrows represent computational steps while boxes
represent the data set. The steps of Weighted Gaussian can run in parallel on
parts of the data set while greedy algorithms must run sequentially over the
data set. Each iteration of SIPS can run in parallel, but the iterations must
be done sequentially.

We benchmark the algorithms in several ways, and we find that DP-SIPS and
Weighted Gaussian scale well with large data sets, while DPSU and GW do not.

For the scalability experiments, we implement all of the algorithms in PySpark
to take advantage of parallelism within the algorithms. We then run the
algorithms on Amazon Elastic Map Reduce (EMR) clusters, and we tailor the
PySpark session settings to the number of cores and memory allocation for both
driver and executor nodes to reflect the available resources of the machines.

To benchmark the algorithms, we measure the amount of time required to run the
algorithm (after the dataset has already been read in to PySpark) and return
the number of partitions that were discovered, in order to ensure PySpark’s
lazy evaluation executed the algorithm during the timing phase. See Appendix C
for information about the machine specifications.

For our first scalability experiment, we run the algorithms on subsets of
Synthetic_80M, a synthetic data set with 80 million users and 4.6 billion
items in total, and benchmark the algorithms as the number of users increases.
Figure 9 shows the results of this experiment. DP-SIPS and Weighted Gaussian
scale well to large data sets while DPSU and GW timeout or run out of memory
on subsets with just 2 million and 10 million users, respectively. We ran this
experiment on an EMR cluster with 8 nodes, each of size m5a.8xlarge (each node
has 32 processor cores and 128 GB of memory).

Figure 9. Algorithm runtimes on subsets of Synthetic_80M of given sizes on a
cluster of 8 m5a.8xlarge nodes. The algorithms were run with the same
hyperparameters as those listed in Table 2. DPSU and GW timed out after
subsets with 2 million and 10 million users, respectively.

To further investigate the scaling behaviors of the algorithms, we ran more
scalability experiments on a smaller data set, Amazon, which consists of
product reviews from 4 million users. Figure 10 shows the results from the
first experiment, in which we measure the runtimes of the algorithms on
clusters with a single master node and varying numbers of core nodes. All
nodes are of type m5a.2xlarge, which has 8 processor cores and 32 GB of
memory. Figure 10 illustrates the following runtime trends:

  * •

The runtimes of Weighted Gaussian and SIPS both decrease as the number of core
nodes increases.

  * •

The runtimes of DPSU and GW remain approximately the same even as the number
of core nodes increases.

These results confirm the intuition that Weighted Gaussian and SIPS are both
parallelizable and thus scale well with increased cluster sizes, even on large
data sets. Additionally, it confirms our observation that since DPSU and GW
need to iterate sequentially over each user, increasing the cluster size does
little to improve their running times.

Figure 10. Algorithm runtimes on Amazon dataset with increasing number of
m5a.2xlarge cores in the cluster. The algorithms were run with the same
hyperparameters as those listed in Table 2.

We run an additional scalability experiment on Amazon: instead of increasing
the number of core nodes, we increase the sizes of the core nodes. We use a
single master node and two core nodes, and increase the sizes of all three
nodes. Table 11 shows similar trends as the previous experiment: Weighted
Gauss and SIPS scale with increased node sizes while DPSU and GW do not. So,
even when the sizes of the machines increase, this makes little difference in
the runtimes of DPSU and SIPS.

Figure 11. Algorithm runtimes on Amazon dataset with increasing sizes of nodes
in a cluster with 1 master node and 2 core nodes. The algorithms were run with
the same hyperparameters as those listed in Table 2.

Comparing Figures 10 and 11, we see that for DP-SIPS, the communication
overhead between machines is small. Specifically we can compare SIPS’s
performance on 8 cores in Figure 10, a cluster with 1 master node and 8 core
nodes each of type m5a.2xlarge (8 CPU cores, 32 GB memory) to m5a.8xlarge (32
CPU cores, 128 GB memory) in Figure 11, a cluster with 1 master node and 2
core nodes. In both, the core nodes have in total $8\cdot 8=2\cdot 32$ CPU
cores and $8\cdot 32=2\cdot 128$ GB of memory among them. The runtimes of SIPS
in the first case is $58.41$ seconds while in the second case, it is $66.28$
seconds, which suggests that the communication overhead is small. In fact,
SIPS runs faster on a cluster of 9 smaller machines than a cluster of 3 larger
machines. This discrepancy could be the result of any number of factors in the
machine specifications.

Our takeaway is clear: the scalability experiments confirm that DPSU and GW
are not suitable for massive data sets and the prohibitively poor run times of
DPSU and GW on industrial-scale data sets cannot be overcome by increasing the
computational power of the cluster (either in number of machines or even
machine size) because the algorithms are not parallelizable on a fundamental
level. For larger data sets, SIPS and Weighted Gaussian are the only feasible
algorithms and SIPS consistently has improved accuracy over Weighted Gaussian.

##  6\. Discussion

Differentially private partition selection (or set union) is a fundamental
problem for many private data analysis tasks. Prior approaches to this problem
either suffer from poor accuracy or are prohibitively slow on large data sets
because they are designed to sequentially iterate over the users. We present a
simple algorithm for differentially private partition selection that achieves
accuracy that is comparable to DPSU and runtimes that scale well with
increased computational resources.

###  6.1. Unsuccessful Attempts

One natural question we explored is: can we boost the accuracy of DPSU using
our method of iterating several times and removing partitions from the data
set that were previously returned? We found that empirically this did not
boost the accuracy on the data sets that we tested. Intuitively, this makes
sense because the per-user Gaussian update policy in DPSU prevents users from
adding weight to items that have already reached the buffered threshold. Such
items will very likely be released after the noise addition and thresholding
steps. The DPSU approach achieves the same goal as iterating, albeit with more
precision due to the greedy nature of the Policy Gaussian per-user updates.

###  6.2. Future work

This work raises several natural questions. From a practical perspective, one
may wonder whether there is a scalable algorithm with higher accuracy than DP-
SIPS. Another interesting line of inquiry could consider theoretical
guarantees on the number of partitions released, possibly under some
distributional assumptions on the data. We believe these lines of inquiry
would give important insights into a fundamental problem in private data
analysis.

##  7\. Acknowledgments

We are thankful to the anonymous reviewers for the helpful feedback. This work
was supported by Tumult Labs.

## References

  * [1] Kareem Amin, Jennifer Gillenwater, Matthew Joseph, Alex Kulesza, and Sergei Vassilvitskii.  Plume: Differential privacy at scale.  arXiv preprint arXiv:2201.11603, 2022.
  * [2] Skye Berghel, Philip Bohannon, Damien Desfontaines, Charles Estes, Sam Haney, Luke Hartman, Michael Hay, Ashwin Machanavajjhala, Tom Magerlein, Gerome Miklau, Amritha Pai, William Sexton, and Ruchit Shrestha.  Tumult Analytics: a robust, easy-to-use, scalable, and expressive framework for differential privacy.  arXiv preprint arXiv:2212.04133, December 2022.
  * [3] Bot_Developer.  Daily financial news for 6000+ stocks.  https://www.kaggle.com/datasets/miguelaenlle/massive-stock-news-analysis-db-for-nlpbacktests.
  * [4] Mark Bun and Thomas Steinke.  Concentrated differential privacy: Simplifications, extensions, and lower bounds.  In Theory of Cryptography Conference, pages 635–658. Springer, 2016\.
  * [5] Clément L Canonne, Gautam Kamath, and Thomas Steinke.  The Discrete Gaussian for Differential Privacy.  In Advances in Neural Information Processing Systems, volume 33, pages 15676–15688. Curran Associates, Inc., 2020.
  * [6] Ricardo Silva Carvalho, Ke Wang, and Lovedeep Singh Gondara.  Incorporating item frequency for differentially private set union.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 36, pages 9504–9511, 2022.
  * [7] Cynthia Dwork, Krishnaram Kenthapadi, Frank McSherry, Ilya Mironov, and Moni Naor.  Our data, ourselves: Privacy via distributed noise generation.  In Annual international conference on the theory and applications of cryptographic techniques, pages 486–503. Springer, 2006.
  * [8] Cynthia Dwork, Frank McSherry, Kobbi Nissim, and Adam Smith.  Calibrating noise to sensitivity in private data analysis.  In Theory of cryptography conference, pages 265–284. Springer, 2006\.
  * [9] Google.  Privacy on Beam.  https://github.com/google/differential-privacy/tree/main/privacy-on-beam, November 2022.
  * [10] Google and OpenMined.  PipelineDP.  https://pipelinedp.io, November 2022.
  * [11] Sivakanth Gopi, Pankaj Gulhane, Janardhan Kulkarni, Judy Hanwen Shen, Milad Shokouhi, and Sergey Yekhanin.  Differentially private set union.  In International Conference on Machine Learning, pages 3627–3636. PMLR, 2020.
  * [12] Kritanjali Jain.  Amazon reviews.  https://www.kaggle.com/datasets/ kritanjalijain/amazon-reviews.
  * [13] Gautam Kamath, Jerry Li, Vikrant Singhal, and Jonathan Ullman.  Privately learning high-dimensional distributions.  In Conference on Learning Theory, pages 1853–1902. PMLR, 2019.
  * [14] Aleksandra Korolova, Krishnaram Kenthapadi, Nina Mishra, and Alexandros Ntoulas.  Releasing search queries and clicks privately.  In Proceedings of the 18th international conference on World wide web, pages 171–180, 2009.
  * [15] Tumult Labs.  Tumult Analytics.  https://tmlt.dev, December 2022.
  * [16] Sebastian Meiser and Esfandiar Mohammadi.  Tight on budget? tight bounds for r-fold approximate differential privacy.  In Proceedings of the 2018 ACM SIGSAC Conference on Computer and Communications Security, pages 247–264, 2018.
  * [17] Lakshmipathi N.  Imdb dataset of 50k movie reviews.  https://www.kaggle.com/datasets/lakshmi25npathi/imdb-dataset-of-50k-movie-reviews.
  * [18] Judy Hanwen Shen.  Ask reddit.  https://github.com/heyyjudes/differentially-private-set-union/tree/ea7b39285dace35cc9e9029692802759f3e1c8e8/data.
  * [19] Thought Vector and Stuart Axelbrooke.  Customer support on twitter.  https://www.kaggle.com/datasets/thoughtvector/customer-support-on-twitter.
  * [20] Mark Wijkhuizen.  Simple/normal wikipedia abstracts v1.  https://www.kaggle.com/datasets/markwijkhuizen/simplenormal-wikipedia-abstracts-v1.
  * [21] Royce J Wilson, Celia Yuxin Zhang, William Lam, Damien Desfontaines, Daniel Simmons-Marengo, and Bryant Gipson.  Differentially private SQL with bounded user contribution.  Proceedings on Privacy Enhancing Technologies, 2:230–250, 2020\.

##  Appendix A Proof of Theorem 1

###### Proof of Theorem 1.

Since the algorithm is simply computing a stable histogram, the proof follows
from standard arguments.

Let us denote Algorithm 2 by $\mathcal{M}$. Fix any $\rho>0$, any
$\delta\in(0,1)$, and any $\Delta_{0}\in\mathbb{N}$. Fix two neighboring data
sets $x=(W_{1},\ldots,W_{N})$ and $x^{\prime}$, where $x^{\prime}$ contains
one extra user set $W^{*}$. Let $E$ be the event that
$\mathcal{M}(x^{\prime})\subseteq\cup_{i\in[N]}\overline{W_{i}}$; that is, the
partitions returned for data set $x^{\prime}$ are in the support of
$\mathcal{M}(x)$. We will first argue that, conditioned on event $E$, the
output distributions of $\mathcal{M}(x)$ and $\mathcal{M}(x^{\prime})$ satisfy
pure zCDP (Definition 4). Then we will argue that, because of the thresholding
step, event $E$ occurs with probability at least $1-\delta$.

If we condition both $\mathcal{M}(x)$ and $\mathcal{M}(x^{\prime})$ on event
$E$, then by the Gaussian mechanism (Definition 12) and post-processing (Lemma
6), Algorithm 2 satisfies pure zCDP (Definition 4).

Now, it remains to show that $\Pr[E]\geq 1-\delta$. First, note that each user
contributes to most $\Delta_{0}$ items in the histogram, and further note that
event $E^{c}$ occurs when at least one item from $W^{*}$ that is not in $x$ is
released. Let $W^{\prime}$ denote the set of items in $W^{*}$ that do not
appear in $x$. Then,

$\displaystyle\Pr[E]$ $\displaystyle=\Pr[\forall u\in
W^{\prime},u\notin\mathcal{M}(x)]$ $\displaystyle=\Pr[\cap_{u\in
W^{\prime}}\hat{H}[u]\leq T]$ $\displaystyle=\Pr[\cap_{u\in
W^{\prime}}H[u]+Z_{u}\leq T]\quad\text{For $Z_{u}\sim\mathcal{N}(0,1/2\rho)$}$
$\displaystyle=\prod_{u\in W^{\prime}}\Pr[H[u]+Z_{u}\leq T]\quad\text{By
independence of $Z_{u}$'s}$ $\displaystyle=\prod_{u\in
W^{\prime}}\Pr\left[\frac{1}{\sqrt{|W^{\prime}|}}+Z_{u}\leq T\right]$
$\displaystyle=\left(\Pr\left[\frac{1}{\sqrt{|W^{\prime}|}}+Z_{u}\leq
T\right]\right)^{|W^{\prime}|}\quad\text{Since $Z_{u}$'s are i.i.d.}$
$\displaystyle\geq\min_{k\in[\Delta_{0}]}\left\\{\left(\Pr\left[\frac{1}{\sqrt{k}}+Z_{u}\leq
T\right]\right)^{k}\right\\}\quad\text{Since $|W^{\prime}|\leq\Delta_{0}$}$
$\displaystyle\geq 1-\delta\quad\text{By definition of $T$}.$

Therefore, $\mathcal{M}$ is $\delta$-approximate $\rho$-zCDP. ∎

##  Appendix B Privacy Parameter Conversions

Tables 4 and 5 give the conversions between zCDP and DP that we use for the
experiments, as well as the values of $\alpha$ that are used in Corollary 9 to
do the conversions.

$\rho$ | $\delta_{CDP}$ | $\varepsilon$ | $\delta_{DP}$ | $\alpha$
---|---|---|---|---
0.001 | $1\times 10^{-5}$ | 0.14 | $5.00\times 10^{-5}$ | 77.033
0.005 | $1\times 10^{-5}$ | 0.338 | $5.08\times 10^{-5}$ | 37.037
0.01 | $1\times 10^{-5}$ | 0.495 | $4.99\times 10^{-5}$ | 27.128
0.05 | $1\times 10^{-5}$ | 1.2 | $4.99\times 10^{-5}$ | 13.283
0.1 | $1\times 10^{-5}$ | 1.765 | $4.96\times 10^{-5}$ | 9.86
0.5 | $1\times 10^{-5}$ | 4.41 | $4.90\times 10^{-5}$ | 5.127
Table 4. Values of $\rho,\delta_{CDP},\varepsilon,$ and $\delta_{DP}$ used for experiments with fixed $\delta_{DP}$ and varying $\varepsilon$. $\rho$ | $\delta_{CDP}$ | $\varepsilon$ | $\delta_{DP}$ | $\alpha$
---|---|---|---|---
0.005 | $1\times 10^{-9}$ | 0.62 | $1.04\times 10^{-9}$ | 64.073
0.0055 | $1\times 10^{-8}$ | 0.62 | $1.02\times 10^{-8}$ | 58.443
0.006 | $1\times 10^{-7}$ | 0.62 | $1.01\times 10^{-7}$ | 53.732
0.007 | $1\times 10^{-6}$ | 0.62 | $1.01\times 10^{-6}$ | 46.334
0.0083 | $1\times 10^{-5}$ | 0.62 | $1.01\times 10^{-5}$ | 39.398
0.01 | $1\times 10^{-4}$ | 0.62 | $1.01\times 10^{-4}$ | 33.037
0.013 | $1\times 10^{-3}$ | 0.62 | $1.01\times 10^{-3}$ | 25.863
Table 5. Values of $\rho,\delta_{CDP},\varepsilon,$ and $\delta_{DP}$ used for
experiments with fixed $\varepsilon$ and varying $\delta_{DP}$.

##  Appendix C Amazon EMR Cluster Specifications

Table 6 shows the number of CPU cores and amount of memory available to each
type of machine on Amazon EMR that we use in our experiments.

Machine Size | Cores | Memory (GB)
---|---|---
m5a.xlarge | 4 | 16
m5a.2xlarge | 8 | 32
m5a.4xlarge | 16 | 64
m5a.8xlarge | 32 | 128
Table 6. Number of cores and memory per EMR machine of each size.


# A Non-Negative Matrix Factorization Game

Satpreet H. Singh

###### Abstract

We present a novel game-theoretic formulation of Non-Negative Matrix
Factorization (NNMF), a popular data-analysis method with many scientific and
engineering applications. The game-theoretic formulation is shown to have
favorable scaling and parallelization properties, while retaining
reconstruction and convergence performance comparable to the traditional
Multiplicative Updates (Lee & Seung, 1999) algorithm.

Game Theory, Matrix Factorization, Machine Learning, Data Science

##  1 Introduction

Non-Negative Matrix Factorization (NNMF) is a matrix decomposition method that
approximates a low-rank non-negative matrix $\mathbf{X}\in\mathbb{R}_{\geq
0}^{I\times J}$, using two non-negative matrices
$\mathbf{W}\in\mathbb{R}_{\geq 0}^{I\times K},$ and
$\mathbf{H}\in\mathbb{R}_{\geq 0}^{K\times J},$ of predetermined rank $K$
according to:

$\displaystyle\mathbf{X}\simeq\mathbf{W}\mathbf{H}$ (1)

The non-negative rank $K$ of the non-negative matrix $\mathbf{X}$ can be
higher than the traditional matrix rank over the real field.

NNMF has a long history of use in both scientific and engineering settings due
to it’s tendency to reveal interesting properties about the underlying data.
For example, Lee & Seung (1999) popularized NNMF after applying it to natural
images, showing that it can learn meaningful visual representations ”the parts
of objects”. NNMF has also been applied to large-scale textual corpora to
learn topic-models (Lee & Seung, 1999; Pauca et al., 2004) and to high-
throughput omics data for time-course analysis (Stein-O’Brien et al., 2018).

Recently Gemp et al. (2020, 2021), reformulated a data-analysis algorithm
called Principal Components Analysis (PCA) (Jolliffe, 2002) as a $K$-player
game, to show that it can be be massively scaled and implemented in a
distributed manner. Inspired by their work, we present a novel game-theoretic
formulation of Non-Negative Matrix Factorization (NNMF) and provide an
empirical analysis of our formulation against a traditional NNMF algorithm
(Lee & Seung, 1999).

##  2 PCA through Nash Equilibrium

Principal Components Analysis (PCA) is a data-analysis algorithm that uses an
orthogonal linear transformation to transform a set of observations
$\mathbf{X}\in\mathbb{R}^{n\times d}$ to a new coordinate system, such that
the projections of the data onto the new coordinates have decreasing variance.
In other words, the first coordinate (known as the first principal component)
is the direction of maximal variance of the data, the second (orthogonal)
coordinate is the direction of second-highest variance of the data, and so on.
[See Jolliffe (2002) for a comprehensive review.]

PCAs are typically calculated in one of two ways: (1) using the eigen-
decomposition of the covariance matrix $\Sigma$ of the data matrix
$\mathbf{X}$ or (2) using the singular value decomposition (SVD) of the data
matrix $\mathbf{X}$.

Using the former method, let the data matrix be
$\mathbf{X}\in\mathbb{R}^{n\times d}$, and it’s covariance matrix be
$M^{d\times d}=\frac{1}{n}\mathbf{X}^{T}\mathbf{X}$. PCA then reduces to the
eigenvalue problem $MV=V\Lambda$, where $\Lambda$ is a diagonal matrix
comprised of eigenvalues (known as PC loadings) and the columns of the
orthonormal matrix $V$ provides the sought PC directions.

###  2.1 Game-theoretic PCA formulation

In a recent work, Gemp et al. (2020) showed how the top-$k$ right singular
vectors of an $n$-observation data matrix $\mathbf{X}\in\mathbb{R}^{n\times
d}$ can be calculated by a $k$-player game. Alternatively, if one is working
with a mini-batch $\mathbf{X_{t}}$ of $n^{\prime}<n$ observations, the top-$k$
right singular vectors are given by the top-$k$ eigenvectors of the positive
semidefinite sample covariance matrix
$\Sigma=\mathbb{E}[\,\frac{1}{n^{\prime}}\mathbf{X_{t}}^{T}\mathbf{X_{t}}]\,$.

In their game formulation, each of $k$ players receives access to a
minibatches of the data and ‘owns’ one of the $k$ (approximate) eigenvectors
$\hat{v}_{i}$. Each player then tries to maximize a utility function
$u_{i}(\hat{v}_{i}|\hat{v}_{j<i})$ that is unique to them:

$\displaystyle\max_{\hat{v}_{i}^{\top}\hat{v}_{i}=1}\Big{\\{}u_{i}(\hat{v}_{i}|\hat{v}_{j<i})$
$\displaystyle=\overbrace{\hat{v}_{i}^{\top}\Sigma\hat{v}_{i}}^{\text{Variance}}-\sum_{j<i}\overbrace{\frac{\langle\hat{v}_{i},\Sigma\hat{v}_{j}\rangle^{2}}{\langle\hat{v}_{j},\Sigma\hat{v}_{j}\rangle}}^{\perp\text{-penalty}}\Big{\\}}$
(2)

Here, all expressions involving the inequality $(j<i)$ indicate that players
with a lower index (‘parents’) have already played before the current
(‘child’) player. As shown above, utility functions consist of two terms that
balance competing objectives: one that tries to maximize the variance captured
by the player’s eigenvector, and the other that tries to maximally
orthogonalize this eigenvector to all parent eigenvectors. To reduce
underspecification, eigenvectors are restricted to lie on the unit sphere
($\hat{v}_{i}^{(t)}\in\mathcal{S}^{d-1}$).

1: Given: data minibatches $\mathbf{X_{t}}\in\mathbb{R}^{m\times d}$, total
iterations $T$, initial vector $\hat{v}_{i}^{0}\in\mathcal{S}^{d-1}$ (unit
sphere in $d-1$ dimensions), and step size $\alpha$.

2: $\hat{v}_{i}\leftarrow\hat{v}_{i}^{0}$

3: for $t=1:T$ do

4: $\nabla_{\hat{v}_{i}}\leftarrow
2\mathbf{X_{t}}^{\top}\Big{[}\mathbf{X_{t}}\hat{v}_{i}-\sum_{j<i}\frac{\langle\mathbf{X_{t}}\hat{v}_{i},\mathbf{X_{t}}\hat{v}_{j}\rangle}{\langle\mathbf{X_{t}}\hat{v}_{j},\mathbf{X_{t}}\hat{v}_{j}\rangle}\mathbf{X_{t}}\hat{v}_{j}\Big{]}$

5:
$\nabla^{R}_{\hat{v}_{i}}\leftarrow\nabla_{\hat{v}_{i}}-\langle\nabla_{\hat{v}_{i}},\hat{v}_{i}\rangle\hat{v}_{i}$

6: $\hat{v}_{i}^{\prime}\leftarrow\hat{v}_{i}+\alpha\nabla^{R}_{\hat{v}_{i}}$

7: $\hat{v}_{i}\leftarrow{\hat{v}_{i}^{\prime}}/{||\hat{v}_{i}^{\prime}||}$

8: broadcast($\hat{v}_{i}$) to all other players

9: end for

10: return $\hat{v}_{i}$

Algorithm 1 PCA game (Gemp et al., 2020)

The complete algorithm is given in Algorithm 1. At each iteration, each player
updates their eigenvector $\hat{v}_{i}$ using a fixed step size Riemannian
gradient ascent update $\nabla^{R}_{\hat{v}_{i}}$, that ensures that each
iterate lies on the acceptable manifold (Bonnabel, 2013). The different
utility functions and order of play results in player-$i$ owning the
eigenvector $\hat{v}_{i}$ associated with the $i^{th}$ largest eigenvalue
$\lambda_{i}$. In addition to the above formulation that assumes ordered play,
Gemp et al. (2020) empirically show that simultaneous play also converges to
Nash equilibrium.

##  3 NNMF as a multiplayer game

Before we formulate NNMF as a game, we describe traditional algorithms used
for it.

###  3.1 Traditional algorithms for NNMF

In it’s most general form, NNMF is an optimization problem that involves
minimizing some loss function $L(\mathbf{X},\mathbf{WH})$ that penalizes the
divergence between the original data matrix $\mathbf{X}$ and it’s
reconstruction $\mathbf{WH}$.

We will work with the following square Frobenius norm penalty for the rest of
this manuscript, though more sophisticated objectives, including those that
include various regularization terms, are possible (Berry et al., 2007):

$\displaystyle\underset{w_{i,k},h_{k,j}\geq
0}{\operatorname{minimize}}||\mathbf{X}-\mathbf{W}\mathbf{H}||^{2}_{F}$ (3)

This optimization formulation is non-convex and obtaining globally optimum
solutions is known to be NP-Hard (Vavasis, 2010). However, interest in NNMF
has led to the development of many empirically successful algorithms for
finding local optima. These algorithms fall into three general classes (see
Berry et al. (2007) for a thorough treatment):

#### 1\. Non-negative Alternating Least squares (NALS):

While the optimization in equation 3.1 is not convex in both $\mathbf{W}$ and
$\mathbf{H}$, it is indeed convex in each term taken individually (biconvex).
NALS solves a least-squares problem alternatively for $\mathbf{W}$ and
$\mathbf{H}$, projecting the iterates into the positive orthant at each
iteration, till convergence is observed.

#### 2\. Projected Gradient (PG):

Similar to NALS, PG performs a gradient ascent step on $\mathbf{W}$ and
$\mathbf{H}$ alternatively, and projects the resulting iterates into the
positive orthant at each iteration, till convergence is observed.

1: Given: data matrix, $\mathbf{X}\in\mathbb{R}^{I\times J}$, latent dimension
$K$, maximum iterations $T$, randomly initialized
$\mathbf{W}^{0}\in\mathbb{R}_{++}^{I\times K}$ and
$\mathbf{H}^{0}\in\mathbb{R}_{++}^{K\times J}$.

2: for $t=1:T$ do

3:
${\displaystyle\mathbf{H}^{t+1}\leftarrow\mathbf{H}^{t}\odot{\frac{((\mathbf{W}^{t})^{T}\mathbf{V})}{((\mathbf{W}^{t})^{T}\mathbf{W}^{t}\mathbf{H}^{t})}}}$

4:
${\displaystyle\mathbf{W}^{t+1}\leftarrow\mathbf{W}^{t}\odot{\frac{(\mathbf{V}(\mathbf{H}^{t+1})^{T})}{(\mathbf{W}^{t}\mathbf{H}^{t+1}(\mathbf{H}^{t+1})^{T})}}}$

5: end for

6: return $\mathbf{W}^{t}$ and $\mathbf{H}^{t}$ Note: $\odot$ and
$\mathbf{X}/\mathbf{Y}$ are elementwise multiplication and division
respectively.

Algorithm 2 Multiplicative updates (Lee & Seung, 1999)

#### 3\. Multiplicative updates (MU):

In PG, the choice of gradient step size is typically done heuristically. As an
improvement, Lee & Seung (1999) popularized a variant of PG that uses an
adaptive step size and results in a multiplicative update for $\mathbf{W}$ and
$\mathbf{H}$ at each iteration. The complete algorithm is given in Algorithm
2.

###  3.2 Game-theoretic NNMF formulation

Figure 1: NNMF as a multiplayer game: NNMF is formulated as a game between $I$
row-players, each owning one row of the $\mathbf{W}$ matrix, that play with
$J$ column-players, each owning one column of the $\mathbf{H}$ matrix. In each
interaction, one row-player $i$ and one column-player $j$ maximize their
respective utilities to better reconstruct one element of the data-matrix
$\mathbf{X}_{[i,j]}$. Additional games between the column players can be
introduced to nudge the game towards desirable solutions, such as orthogonal
rows for $\mathbf{H}$.

To set up NNMF as a game, we define two sets of players, $I$ row-players and
$J$ column-players and their respective utilities. Each of the $I$ row-players
owns a row $\mathbf{W}_{[i,:]}$ of the $\mathbf{W}$ matrix, and each the $J$
column-players owns a column $\mathbf{H}_{[:,j]}$ of the $\mathbf{H}$ matrix.
Therefore, each player owns a length-$K$ non-negative vector.

Taking inspiration from the game formulation of PCA, we now define a utility
function for each player, such that a player’s local actions (i.e. best
response) achieve a global optimization objective. The primary interaction in
this game is between row-players playing column-players and vice-versa. The
utility function in these interactions is given by:

$\displaystyle
u_{(i,j)}=-\mathbf{\ell}(i,j)=-\frac{1}{2}[\,\mathbf{X}_{[i,j]}-(\mathbf{W}_{[i,:]})^{T}\mathbf{H}_{[:,j]}]\,^{2}$
(4)

In other words, the $i^{th}$ row-player and $j^{th}$ column-player cooperate
to try to reconstruct the $(i,j)^{th}$ element of the data matrix
$\mathbf{X}$. The NNMF game algorithm is summarized in Algorithm 3.

The formulation described up to this point can be regarded as a graphical game
on a bipartite graph with row-players and column players comprising the two
disjoint sets. We now relax the bipartite structure to a more general
graphical game where there are interactions between some identical players,
similar to replicator or winner-take-all dynamics from evolutionary game
theory (Sandholm, 2010). Column players (owning columns of the $\mathbf{H}$
matrix) now have self-games as a way to induce certain forms of regularization
that aid interpretation of the resulting matrix $\mathbf{H}$. For example, to
encourage non-overlapping rows in $\mathbf{H}$, we could have a column player
play against itself in such a way that the smallest element of the column-
vector it owns is set to zero. A more aggressive regularization strategy
towards the same goal could be that all but the largest element of the column-
vector stay non-zero after the self-interaction. Henceforth, we call these
strategies J-min and J-max respectively, In practice by instead of setting an
element of the column-vector to zero, we multiply it by $0.99$ to nudge it in
the right direction.

1: Given: data matrix, $\mathbf{X}\in\mathbb{R}^{I\times J}$, latent dimension
$K$, maximum iterations $T$, step size $\eta$, randomly initialized
$\mathbf{W}^{0}\in\mathbb{R}_{++}^{I\times K}$ and
$\mathbf{H}^{0}\in\mathbb{R}_{++}^{K\times J}$.

2: for $t=1:T$ do

3: for $i=1:I$ do

4: for $j=1:J$ do

5:
$\nabla^{t}_{W}\leftarrow-[\,\mathbf{X}_{[i,j]}-(\mathbf{W}^{t}_{[i,:]})^{T}\mathbf{H}^{t}_{[:,j]}]\,(\mathbf{H}^{t}_{[:,j]})^{T}$

6:
$\nabla^{t}_{H}\leftarrow-[\,\mathbf{X}_{[i,j]}-(\mathbf{W}^{t}_{[i,:]})^{T}\mathbf{H}^{t}_{[:,j]}]\,(\mathbf{W}^{t}_{[i,:]})^{T}$

7: $\mathbf{W}^{t+1}_{[i,:]}\leftarrow
max[\,\mathbf{0},\mathbf{W}^{t}_{[i,:]}-\eta\nabla^{t}_{W}]\,$

8: $\mathbf{H}^{t+1}_{[:,j]}\leftarrow
max[\,\mathbf{0},\mathbf{H}^{t}_{[:,j]}-\eta\nabla^{t}_{H}]\,$

9: end for

10: end for

11: end for

12: return $\mathbf{W}^{t}$ and $\mathbf{H}^{t}$

Algorithm 3 NNMF game (simultaneous play)

###  3.3 Parallelization and scaling properties

The game-theoretic formulation of NNMF has several favorable parallelization
and scaling properties:

#### Data sharding:

Each row-player and each column-player only require access to their respective
row or column of the data matrix $\mathbf{X}$. In Algorithm 3 only a single
element of the data matrix is required for every row-column player
interaction. If each row-player and each column-player make a local copy of
their associated row or column of the data-matrix, that would be sufficient
for all possible row-column interactions they have. Furthermore, if we assume
ordered play in row-column interactions, only the leader need hold the data
matrix. They can then share the required ground-truth matrix element with the
other player during the interaction.

In our formulation, we assume that each (row or column) player is on a
separate computational core. However, the algorithm is indifferent to the
actual distribution of players across computational cores. This distribution
can be arbitrary or be dictated by some optimal trade-off between inter-core
communication cost and per-core memory cost.

#### Computational cost:

The computational cost, in floating point operations (FLOPS), of multiplying
one $I\times K$ matrix with a $K\times J$ matrix is $O(IJK)$.

Therefore one MU iteration (Algorithm 2) costs  $O(IJK)$
[${((\mathbf{W}^{t})^{T}\mathbf{V})}$] \+ $O(IK^{2}+JK^{2})$
[${((\mathbf{W}^{t})^{T}\mathbf{W}^{t}\mathbf{H}^{t})}$] + $O(KJ)$
[$\mathbf{X/Y}$] + $O(KJ)$ [$\odot$] for the $\mathbf{H}$ update. A similar
calculation applies to the $\mathbf{W}$ update.

One NNMF iteration (Algorithm 3) involves $IJ$ row-column player interactions,
each with a cost $O(K^{2})$
[$\nabla^{t}_{W}\leftarrow-[\,\mathbf{X}_{[i,j]}-(\mathbf{W}^{t}_{[i,:]})^{T}\mathbf{H}^{t}_{[:,j]}]\,(\mathbf{H}^{t}_{[:,j]})^{T}$]
\+ $O(K)$ [
$\mathbf{W}^{t+1}_{[i,:]}\leftarrow\mathbf{W}^{t}_{[i,:]}-\eta\nabla^{t}_{W}$
] for the $W_{[i,:]}$ update. A similar calculation applies for the
$\mathbf{H}$ update. This is comparable with the cost of the MU update for
small $K$, while being highly parallelizable.

##  4 Empirical analysis

To validate our new game-theoretic formulation of NNMF, we empirically compare
it with the traditional Multiplicative Updates (MU) algorithm using synthetic
datasets.

Figure 2: Synthetic data for experiments: Data matrix
$\mathbf{X_{syn}}\in\mathbb{R}^{I\times J}$ reconstructed from matrix
multiplication of $\mathbf{H_{syn}}\in\mathbb{R}_{++}^{K\times J}$ and
$\mathbf{W_{syn}}\in\mathbb{R}_{++}^{I\times K}$.

#### Synthetic dataset generator:

For all experiments, we set the number and dimension of our observations to be
$I=100$ and $J=20$ respectively, resulting in a data matrix
$\mathbf{X_{syn}}\in\mathbb{R}^{I\times J}$. We set the true number of latent
factors to $K=3$, and generate a random Uniform$(0,1)$ non-negative mixing
matrix $\mathbf{W_{syn}}\in\mathbb{R}_{++}^{I\times K}$. We generate the non-
negative basis matrix $\mathbf{H_{syn}}\in\mathbb{R}_{++}^{K\times J}$ by
applying a Gaussian filter (window-length=3) to the rows of a ${K\times J}$
random matrix with elements drawn from a Uniform$(0,1)$ distribution.
$\mathbf{X_{syn}}$ is then obtained multiplying $\mathbf{W_{syn}}$ and
$\mathbf{H_{syn}}$ (Figure 2).

#### Experiment and Results:

Figure 3: Convergence of algorithm variants [Top]: Game-theoretic variants of
NNMF ($\eta=0.001$) have convergence rates comparable to the traditional
Multiplicative Updates (MU) algorithm.
Ground truth vs. Recovered Bases [Lower]: Basis matrices $\mathbf{H_{syn}}$
recovered from MU and Game variants look similar, up to permutations of order
(colors). The J-min and J-max variants only show slight effects of
regularization, possibly because the ground-truth basis is not exactly
orthogonal.  Figure 4: Reconstruction error of algorithm variants evaluated
over multiple randomly generated synthetic datasets: Game-theoretic NNMF
variants result in reconstruction errors that are not consistently
significantly different compared to reconstruction errors from the traditional
MU algorithm. Error bars represent $\pm$ 1 s. d., calculated over 4 random
initializations. Figure 5: Low dimensional projection of iterate trajectories
for $\mathbf{\hat{H}_{syn}}$ basis matrix rows: Trajectories for
Multiplicative Updates (MU) [orange], NNMF game vanilla [purple], Jmin
[green], and Jmax [blue] variants (using same initialization as in Figure 3).
Each algorithm has $K=3$ trajectories associated with it, one for each row of
the recovered basis matrix $\mathbf{\hat{H}_{syn}}$, each starting at a
$\bullet$ marker and ending at a $\times$ marker. Trajectories are observed to
be relatively smooth for the (adaptive step-size) MU algorithm and relatively
jagged for the (fixed step-size) NNMF game variants. Though iterates seem to
be moving away from each other, it is not clear if this can be regarded as a
signature of congestion as claimed for iterates in the PCA game Gemp et al.
(2020). Low dimensional space is generated from applying PCA on all
$\mathbf{\hat{H}_{syn}}$ basis vectors from all algorithms stacked together.
PC1, PC2 and PC3 explain 69, 15 and 8 percent of the total variance of the
data respectively.

We generate a synthetic dataset and compare its decomposition using MU
algorithm to the one obtained using our game-theoretic NNMF algorithm and it’s
two regularized variants, J-min and J-max, Section 3.2 (using a fixed gradient
step size of $\eta=0.001$ across all NNMF executions). We find that the
convergence trajectories of the MU and NNMF-Game algorithms are comparable,
and that the recovered basis matrices $\mathbf{\hat{H}}_{syn}$ look quite
similar (Figure 3).

Next, we examine the evolution of the non-negative basis (row) vectors that
comprise the basis matrix $\mathbf{\hat{H}}_{syn}$ by plotting a low-
dimensional projection of the trajectories of the iterates (Figure 5). This
reveals trajectories that appear smooth for the (adaptive step size) MU
algorithm and relatively jagged for our (fixed step size) NNMF game variants.

Finally, we compare the reconstruction performance
($Error=||\mathbf{X}-\mathbf{\hat{X}}||_{2}$) of these algorithms across 10
new independently randomly generated synthetic datasets, with 4 randomly
initialized episodes for each. We find that the performance of the game-
theoretic NNMF variants and that of the traditional MU algorithm is not
significantly different (Figure 4).

##  5 Discussion

We present a game formulation of the popular Non-negative Matrix Factorization
algorithm with favorable parallelization and scaling properties. The game
structure admits variants that can induce certain forms of interpretability-
enhancing regularization such as encouraging rows in the recovered basis
matrix $\mathbf{\hat{H}}$ to be orthogonal. We provide empirical evidence that
our algorithms have convergence and reconstruction performance similar to
those of the traditional Multiplicative Updates (MU) (Lee & Seung, 1999)
algorithm. Visual inspection of convergence trajectories for rows of the
recovered basis matrix $\mathbf{\hat{H}}$ reveals that iterates tend to vary
relatively smoothly for the (adaptive step-size) MU algorithm but relatively
jaggedly for our (fixed step-size) NNMF game variants.

For future work, we plan to formally analyse the convergence properties of the
graphical game variants introduced in this paper, and explore adaptive
learning rate and additional regularization strategies that could result in
smoother convergence trajectories and additional desirable properties in the
recovered factor matrices. Exploring extensions to the more general Tensor
Factorization setting, and connections to Message Passing inference algorithms
are other directions that seem promising.

## Acknowledgements

We thank Prof. Lillian Ratliff for helpful discussions.

## References

  * Berry et al. (2007) Berry, M. W., Browne, M., Langville, A. N., Pauca, V. P., and Plemmons, R. J.  Algorithms and applications for approximate nonnegative matrix factorization.  _Computational statistics & data analysis_, 52(1):155–173, 2007.
  * Bonnabel (2013) Bonnabel, S.  Stochastic gradient descent on Riemannian manifolds.  _IEEE Transactions on Automatic Control_ , 58(9):2217–2229, 2013.
  * Gemp et al. (2020) Gemp, I., McWilliams, B., Vernade, C., and Graepel, T.  Eigengame: PCA as a Nash Equilibrium.  _arXiv preprint arXiv:2010.00554_ , 2020.
  * Gemp et al. (2021) Gemp, I., McWilliams, B., Vernade, C., and Graepel, T.  Eigengame unloaded: When playing games is better than optimizing.  _arXiv preprint arXiv:2102.04152_ , 2021.
  * Jolliffe (2002) Jolliffe, I. T.  Principal components in regression analysis.  In _Principal Component Analysis_. Springer, 2002.
  * Lee & Seung (1999) Lee, D. D. and Seung, H. S.  Learning the parts of objects by non-negative matrix factorization.  _Nature_ , 401(6755):788–791, 1999.
  * Pauca et al. (2004) Pauca, V. P., Shahnaz, F., Berry, M. W., and Plemmons, R. J.  Text mining using non-negative matrix factorizations.  In _Proceedings of the 2004 SIAM International Conference on Data Mining_ , pp. 452–456. SIAM, 2004.
  * Sandholm (2010) Sandholm, W. H.  _Population games and evolutionary dynamics_.  MIT press, 2010.
  * Stein-O’Brien et al. (2018) Stein-O’Brien, G. L., Arora, R., Culhane, A. C., Favorov, A. V., Garmire, L. X., Greene, C. S., Goff, L. A., Li, Y., Ngom, A., Ochs, M. F., et al.  Enter the matrix: factorization uncovers knowledge from omics.  _Trends in Genetics_ , 34(10):790–805, 2018.
  * Vavasis (2010) Vavasis, S. A.  On the complexity of nonnegative matrix factorization.  _SIAM Journal on Optimization_ , 20(3):1364–1377, 2010.