# Spatial-spectral mapping to prepare the frequency entangled qudits

Zi-Xiang Yang  Hubei Key Laboratory of Optical Information and Pattern
Recognition, Wuhan Institute of Technology, Wuhan 430205, China  Zi-Qi Zeng
Hubei Key Laboratory of Optical Information and Pattern Recognition, Wuhan
Institute of Technology, Wuhan 430205, China  Ying Tian  Hubei Key Laboratory
of Optical Information and Pattern Recognition, Wuhan Institute of Technology,
Wuhan 430205, China  Shun Wang  Hubei Key Laboratory of Optical Information
and Pattern Recognition, Wuhan Institute of Technology, Wuhan 430205, China
Ryosuke Shimizu  The University of Electro-Communications, 1-5-1 Chofugaoka,
Chofu, Tokyo, Japan  Hao-Yu Wu<EMAIL_ADDRESS>Hubei Key Laboratory of
Optical Information and Pattern Recognition, Wuhan Institute of Technology,
Wuhan 430205, China  Shilong Liu<EMAIL_ADDRESS>Department of
Physics, University of Ottawa, 25 Templeton Street, K1N 6N5, Ottawa, ON,
Canada  femtoQ Lab, Engineering Physics Department, Polytechnique Montréal,
Montréal, Québec H3T 1JK, Canada  Rui-Bo Jin<EMAIL_ADDRESS>Hubei Key
Laboratory of Optical Information and Pattern Recognition, Wuhan Institute of
Technology, Wuhan 430205, China  Guangdong Provincial Key Laboratory of
Quantum Science and Engineering, Southern University of Science and
Technology, Shenzhen 518055, China

###### Abstract

Entangled qudits, the high-dimensional entangled states, play an important
role in the study of quantum information. How to prepare entangled qudits in
an efficient and easy-to-operate manner is still a challenge in quantum
technology. Here, we demonstrate a method to engineer frequency entangled
qudits in a spontaneous parametric downconversion process. The proposal
employs an angle-dependent phase-matching condition in a nonlinear crystal,
which forms a classical-quantum mapping between the spatial (pump) and
spectral (biphotons) degrees of freedom. In particular, the pump profile is
separated into several bins in the spatial domain, and thus shapes the down-
converted biphotons into discrete frequency modes in the joint spectral space.
Our approach provides a feasible and efficient method to prepare a high-
dimensional frequency entangled state. As an experimental demonstration, we
generate a three-dimensional entangled state by using a homemade variable slit
mask.

††preprint: APS/123-QED

##  I Introduction

Quantum entangled states are serving as essential resources in quantum
technologies, e.g., quantum computation Knill2001 , communications
gisin2007quantum , and measurements clerk2010introduction . The high-
dimensional entangled states (or named as entangled qudits with a dimension
number of $d$) demonstrate significant progress in the aforementioned quantum
applications Erhard2020 . In quantum communication, a high-dimensional quantum
state could carry more information, thus increasing the channel capacities and
also the noise resilience Chen2021 ; dada2011experimental . In quantum
computation, entangled qudits not only have a larger state space to store and
process information but also have the ability to do multiple control
operations simultaneously Lu2019 . These features are significant in the
reduction of the circuit complexity and the acceleration of the algorithm
Reimer2018 . In quantum measurement, a strong reduction in the number of
operations can be achieved by using qudit systems satisfying a certain
relation between their dimensionality and topology Giovannetti2004 . In
addition, entangled qudits are crucial for studies of fundamental quantum
mechanics, i.e., the non-locality in high-dimensional quantum systems
dada2011experimental .

The high-dimensional entangled state could be realized in various degrees of
freedom of the photon, including time LagoRivera2021 , frequency Kues2017 ;
Chen2021 ; Imany2018 , hybrid time-frequency modes Reimer2018 ; Xiang2020 ,
paths Lu2020 , and orbital angular momentum (OAM) liu2018coherent . Regarding
the time-frequency entangled qudits, they are intrinsically suitable for long-
distance transmission in optical fibers, waveguides, as well as free space.
Therefore, it has attracted much attention in recent years Chen2021 ;
Lingaraju:19 ; PhysRevLett.123.123603 ; Bernhard2013 ;
XieZhongShresthaEtAl2015 ; KuesReimerRoztockiEtAl2017 ; Jin2016QST ;
Jin2021APLphoton ; Lixiaoying2020APL ; Morrison2022APL ; Yang2020  .

Generally speaking, the previous methods for the generation of the frequency-
entangled qudits could be classified into four categories: (1) using a linear
pulse shaper, which may include two optical gratings and one spatial light
modulator (SLM) Bernhard2013 ; Liu2022SM ; (2) utilizing a nonlinear cavity,
such as a Fabry-Perot cavity or ring-cavity XieZhongShresthaEtAl2015 ;
KuesReimerRoztockiEtAl2017 ; (3) employing an interferometer, for example, the
spectrally resolved Hong-Ou-Mandel(HOM) interferometer Jin2016QST , quantum
optical synthesis in a dual pump interferometer Jin2021APLphoton , or a
nonlinear interferometer in fibers Lixiaoying2020APL ; (4) developing quantum
state engineering in a nonlinear material, such as a customized poling
nonlinear crystal Morrison2022APL . All the above schemes require
sophisticated modulation devices or specially designed crystals, which may
lead to a large loss or high cost during the preparation. Therefore, it’s
meaningful to explore one simple and efficient method to generate a high-
dimensional frequency entangled state.

Figure 1: The principle of spatial-spectral mapping. (a-c) The different
incident angles of the pump beam correspond to discrete spectral modes. (d)
The concept of combining three beams to prepare a three-dimensional entangled
state. (e) The configuration of the incident angle $\theta$: the angle between
the pump laser and the optical axis of the crystal. (f) The incident angle
distribution on the cross section of the pump beam. (g) The intensity
distribution on the cross section of the pump beam. Figure 2: Simulation of
the spatial-spectral mapping to prepare high-dimensional entangled states.
(a1-a5): The cross section of the pump beam, spatially filtered by different
masks. (b1-b5): The calculated JSIs under the corresponding spatial masks.
(c1-c5): The spectra of the signal photon are obtained by projecting the JSI
onto the vertical axis. See supplemental document for more simulations.

Here, we propose and experimentally demonstrate a feasible method to generate
frequency-entangled qudits using a $\beta$-barium borate (BBO) crystal. This
is possible by employing two improvements in our regime. On the one hand, we
employ an angle-dependent phase-matching (also called birefringent phase-
matching) condition of the BBO crystal, to form a classical-quantum mapping
between the spatial (pump) and spectral (photon pairs) degrees of freedom. On
the other hand, we build a homemade spatial mask, which could make a spatial
modulation on the pump profile. For example, it could separate the pump into
several spatial bins and thus on-demand control the entangled spectral modes
in the joint spectral space. The proposed regime is of great feasibility and
provides one efficient way to prepare the high-dimensional entangled qubit in
the frequency domain.

Figure 3:  (a) The experimental setup and results for preparation of the
frequency entangled qudits. L1(2): lens, LPFs: long-pass filters, PBS:
polarizing beam splitter, M1: mirror, D1(2): detector, TIA: time interval
analyzer, SMFC: single-mode-fiber coupler. (b-g) The measured spectra of the
signal photon using different masks, which are shown in the upper-right corner
of each figure.

##  II Theory and simulation

The biphoton state generated in a spontaneous parametric down conversion
(SPDC) process can be expressed asMosley2008PRL :

$\left|\psi\right\rangle=\int_{0}^{\infty}{\int_{0}^{\infty}{d{\omega_{s}}d{\omega_{i}}f({\omega_{s}},{\omega_{i}},\theta)}}{\hat{a}^{\dagger}}_{s}({\omega_{s}}){\hat{a}^{\dagger}}_{i}({\omega_{i}})\left|0\right\rangle,$
(1)

where the subscripts $s$ and $i$ represent the signal and idler photon,
respectively; $\omega$ is the angular frequency; ${\hat{a}^{\dagger}}$ is the
creation operator; $f({\omega_{s}},{\omega_{i}},\theta)$ is the joint spectral
amplitude (JSA) of the biphoton, where $\theta$ is the incident angle of the
pump, i.e., the angle between the pump laser and the optical axis of the
crystal (See Fig. 1 (e)). $f({\omega_{s}},{\omega_{i}},\theta)$ is the product
of the pump-envelope function $\alpha({\omega_{s},\omega_{i}})$ and the phase-
matching function $\phi({\omega_{s}},{\omega_{i},\theta})$ Mosley2008PRL ,
i.e.,

$f({\omega_{s}},{\omega_{i}},\theta)=\alpha({\omega_{s}},{\omega_{i}})\times\phi({\omega_{s}},{\omega_{i}},\theta).$
(2)

The widely used pump-envelope function is a Gaussian function, i.e.,

$\alpha({\omega_{s}},{\omega_{i}})=\exp\left[-\frac{1}{2}{\left(\frac{{{\omega_{s}}+{\omega_{i}}-{\omega_{p}}}}{{{\sigma_{p}}}}\right)^{2}}\right],$
(3)

where $\omega_{p}$ and $\sigma_{p}$ are the central frequency and bandwidth of
the pump, respectively. The phase-matching function in a nonlinear crystal can
be expressed as Mosley2008PRL

$\phi({\omega_{s}},{\omega_{i}},\theta)={\rm{sinc}}\left(\frac{{\Delta
kL}}{2}\right)\exp\left(\frac{{i\Delta kL}}{2}\right),$ (4)

where $L$ is the crystal length; $\Delta k={k_{s}}+{k_{i}}-{k_{p}}$ is the
wave vector difference between the signal, idler, and pump, respectively. The
wave vector is a function of the incident angle $\theta$.

The principle for generating qudits is shown in Fig. 1. Thanks to the angle
phase-matching (also called birefringent phase-matching) condition, the
different incident angle $\theta_{j}$ corresponds to the different discrete
spectral mode in the joint spectral intensity (JSI)
$|f({\omega_{s}},{\omega_{i}},\theta_{j}){|^{2}}$. For example, Fig. 1 (a-c)
shows the cases of $\theta_{j}=41.29^{\circ}$, $41.79^{\circ}$, and
$42.29^{\circ}$, respectively. The corresponding central wavelengths are 794
nm (827 nm), 810 nm (810 nm), and 827 nm (794 nm) for the signal (idler)
photon, respectively. If we shape the pump profile to $d$ spatial rays, the
output state will become a separate but superposed $d$-dimensional entangled
state:

$\left|\psi\right\rangle=\sum_{j=1}^{d}c_{j}\left|\omega_{s},\omega_{i}\right\rangle_{j}=\sum_{j=1}^{d}c_{j}\left|j\right\rangle_{\omega_{s}}\left|j\right\rangle_{\omega_{i}},$
(5)

Here, $c_{j}$ is the coefficient that depends on both the spatial pump beam
and JSA (JSI). In this case, the JSI could be expressed as:

$JSI=\left|\sum_{j=1}^{d}f({\omega_{s}},{\omega_{i}},\theta_{j})\right|^{2}.$
(6)

As a demonstration, Fig. 1(d) shows three-dimensional entangled qudits.

Under the paraxial approximation, the incident angle of the pump beam is only
determined by the vertical axis. In other words, all the positions on the same
horizontal line have the same incident angle, as shown in Fig. 1(e). This
phenomenon can be intuitively understood in the following way: since the
focusing range (e.g., 50 mm in our scheme) is much larger than the focused
beam radius (less than 1 mm) on the BBO crystal, the contribution from the
horizontal angle can be neglected. Fig. 1(f, g) shows the incident angle and
intensity distribution on the cross section of the pump beam.

Based on the incident angle distribution in Fig. 1(f), we can shape the pump
beam using homemade masks. Fig. 2(a1-a5) and (b1-b5) show the simulations,
where the different masks correspond to the unique JSIs. When no slit impinges
on the pump beam, a tilt and long JSI can be generated, as shown in Fig. 2(a1,
b1). When the center position of the pump beam is blocked using a mask in Fig.
2(a2), the center position in the JSI is lost and a two-mode frequency
entangled state is generated, as shown in Fig. 2(b2). When the two sections
between the center are blocked using the mask in Fig. 1(a3), two related
sections in the JSI are also lost, resulting in a three-dimensional qudit in
Fig. 1(b3). Following a similar way, four- and five-dimensional qudits in Fig.
2(b4, b5) could also be generated using three- and four-line masks in Fig.
2(a4, a5).

By projecting the two-dimensional distribution diagram of JSI onto the
vertical and horizontal axes, the spectral distribution of the signal and
idler photon can be obtained. Fig. 2(c1-c5) shows the corresponding
distribution of the signal photon, which has a distribution of 1, 2, 3, 4, and
5 peaks. The full-width at half-maximum (FWHM) of the peak in Fig. 2(c1) is
19.7 nm. Note by adjusting the width and the center position of the mask,
different JSIs and spectra can be prepared. See the supplemental document for
more simulations. Next, we verify this scheme in an experiment.

##  III Experiment and Results

The experimental setup for generating frequency entangled qudits is shown in
Fig. 3 (a). The laser used in this experiment is a single-transverse-mode and
multi-longitudinal-mode laser diode (LD), which has a central wavelength of
405 nm and a bandwidth of 0.53 nm Cai2022 . The laser beam (with a diameter of
around 2 mm) is spatially filtered by a mask and then focused using a lens
(L1, f=50 mm). The 5-mm-long BBO crystal is designed for a type-II phase-
matched (e$\rightarrow$o+e) SPDC at 810 nm (cut angle: $\theta=41.9^{\circ}$
and $\varphi=0^{\circ}$). The downconverted biphotons are collimated by the
second lens (L2, f=50 mm) and then filtered by a set of long-pass filters
(LPFs). Here, we define the reflected (transmitted) photon from the PBS as the
signal (idler). After being separated by a polarization beam splitter (PBS),
the biphotons are coupled into the single-mode fibers (SMFs), which are
connected to two single-photon detectors D1 and D2 (SPCM-AQRH-10-FC from
Excelitas) and a time interval analyzer (TIA, Picoharp300 from Pico Quanta
Co.). When the pump power is set to 40 mW, we obtain single counts of 762 kHz
for the signal and 325 kHz for the idler. Here, the gap in single counts is
due to the different coupling efficiencies of each channel. The measured
coincidence count is 4.1 kHz.

After the photon counting test, the biphotons’ spectra are measured by a
single-photon level spectrometer (SP2300, Princeton Instrument Co.). Figure
3(b) is the spectrum of the signal photon using a homemade mask with no block
line inside. The FWHM of this spectrum is 54.2 nm. This mask is made of
photocurable resin using a 3D printer. By inserting a mask with one block
line, we obtain the spectrum of the two-dimensional qudits, as shown in Fig.
3(c, d, e, f). To investigate the orthogonality of the frequency qudits, we
set the width of the center block line as 0.3 mm, 0.6 mm, 1.0 mm, and 1.5 mm
in Fig. 3(c, d, e, f), respectively. It can be noticed that the two peaks in
Fig. 3(c, d) are partially overlapped, while the two peaks in Fig. 3(e, f) are
completely separated, indicating that the two frequency components are
orthogonal. The spacing between the two peak centers is 41.1 nm, 49.5 nm, 59.0
nm, and 77.6 nm, respectively. Figure 3(g) depicts the spectrum of the three-
dimensional entangled qudits using a mask with two block lines inside. The
width of each block line is 0.8 mm and the spacing between the two lines is
0.1 mm. The three frequency components in Fig 3(g) are also well separated.

##  IV Conclusion

In conclusion, the present scheme may provide a useful platform to engineer
the more complex frequency entangled state. For example, by employing a
flattop spectrum of the pump, the regime could generate a high-dimensional
maximum entangled state liu2020increasing . By selecting one spatial mask with
a high resolution, i.e., an SLM, the dimensional space of the entangled state
could be further improved. Also, a high-dimensional Bell state may be possible
to perform with the help of the pump modulation technique liu2018coherent . In
addition, the reported spatial-spectral mapping scheme is beneficial to study
the spatial-temporal entangled states Brecht2015PRX .

## Acknowledgments

This work was supported by the National Natural Science Foundations of China
(Grant Numbers 91836102, 12074299, 11704290, and 11904112); Guangdong
Provincial Key Laboratory (Grant No. GKLQSE202102) and Natural Science
Foundation of Hubei Province (2022CFA039).

## References

  * (1) E. Knill, R. Laflamme, and G. J. Milburn, “A scheme for efficient quantum computation with linear optics,” Nature 409, 46–52 (2001).
  * (2) N. Gisin and R. Thew, “Quantum communication,” Nat. Photon. 1, 165–171 (2007).
  * (3) A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt, and R. J. Schoelkopf, “Introduction to quantum noise, measurement, and amplification,” Rev. Mod. Phys. 82, 1155 (2010).
  * (4) M. Erhard, M. Krenn, and A. Zeilinger, “Advances in high-dimensional quantum entanglement,” Nature Reviews Physics 2, 365–381 (2020).
  * (5) Y. Chen, S. Ecker, L. Chen, F. Steinlechner, M. Huber, and R. Ursin, “Temporal distinguishability in Hong-Ou-Mandel interference for harnessing high-dimensional frequency entanglement,” npj Quantum Inform. 7 (2021).
  * (6) A. C. Dada, J. Leach, G. S. Buller, M. J. Padgett, and E. Andersson, “Experimental high-dimensional two-photon entanglement and violations of generalized bell inequalities,” Nature Phys. 7, 677–680 (2011).
  * (7) H.-H. Lu, Z. Hu, M. S. Alshaykh, A. J. Moore, Y. Wang, P. Imany, A. M. Weiner, and S. Kais, “Quantum phase estimation with time-frequency qudits in a single photon,” Advanced Quantum Technologies 3, 1900074 (2019).
  * (8) C. Reimer, S. Sciara, P. Roztocki, M. Islam, L. R. Corté, Y. Zhang, B. Fischer, S. Loranger, R. Kashyap, A. Cino, S. T. Chu, B. E. Little, D. J. Moss, L. Caspani, W. J. Munro, J. Azańa, M. Kues, and R. Morandotti, “High-dimensional one-way quantum processing implemented on d-level cluster states,” Nature Phys. 15, 148–153 (2018).
  * (9) V. Giovannetti, S. Lloyd, and L. Maccone, “Quantum-enhanced measurements: Beating the standard quantum limit,” Science 306, 1330–1336 (2004).
  * (10) D. Lago-Rivera, S. Grandi, J. V. Rakonjac, A. Seri, and H. de Riedmatten, “Telecom-heralded entanglement between multimode solid-state quantum memories,” Nature 594, 37–40 (2021).
  * (11) M. Kues, C. Reimer, P. Roztocki, L. R. Cortés, S. Sciara, B. Wetzel, Y. Zhang, A. Cino, S. T. Chu, B. E. Little, D. J. Moss, L. Caspani, J. Azaña, and R. Morandotti, “On-chip generation of high-dimensional entangled quantum states and their coherent control,” Nature 546, 622–626 (2017).
  * (12) P. Imany, J. A. Jaramillo-Villegas, O. D. Odele, K. Han, D. E. Leaird, J. M. Lukens, P. Lougovski, M. Qi, and A. M. Weiner, “50-GHz-spaced comb of high-dimensional frequency-bin entangled photons from an on-chip silicon nitride microresonator,” Opt. Express 26, 1825–1840 (2018).
  * (13) X. Xiang, R. Dong, R. Quan, Y. Jin, Y. Yang, M. Li, T. Liu, and S. Zhang, “Hybrid frequency-time spectrograph for the spectral measurement of the two-photon state,” Opt. Lett. 45, 2993–2996 (2020).
  * (14) L. Lu, L. Xia, Z. Chen, L. Chen, T. Yu, T. Tao, W. Ma, Y. Pan, X. Cai, Y. Lu, S. Zhu, and X.-S. Ma, “Three-dimensional entanglement on a silicon chip,” npj Quantum Inform. 6 (2020).
  * (15) S. Liu, Z. Zhou, S. Liu, Y. Li, Y. Li, C. Yang, Z. Xu, Z. Liu, G. Guo, and B. Shi, “Coherent manipulation of a three-dimensional maximally entangled state,” Phys. Rev. A 98, 062316 (2018).
  * (16) N. B. Lingaraju, H.-H. Lu, S. Seshadri, P. Imany, D. E. Leaird, J. M. Lukens, and A. M. Weiner, “Quantum frequency combs and Hong-Ou-Mandel interferometry: the role of spectral phase coherence,” Opt. Express 27, 38683–38697 (2019).
  * (17) V. V. Orre, E. A. Goldschmidt, A. Deshpande, A. V. Gorshkov, V. Tamma, M. Hafezi, and S. Mittal, “Interference of temporally distinguishable photons using frequency-resolved detection,” Phys. Rev. Lett. 123, 123603 (2019).
  * (18) C. Bernhard, B. Bessire, T. Feurer, and A. Stefanov, “Shaping frequency-entangled qudits,” Phys. Rev. A 88, 032322 (2013).
  * (19) Z. Xie, T. Zhong, S. Shrestha, X. Xu, J. Liang, Y.-X. Gong, J. C. Bienfang, A. Restelli, J. H. Shapiro, F. N. C. Wong, and C. Wei Wong, “Harnessing high-dimensional hyperentanglement through a biphoton frequency comb,” Nat. Photon. 9, 536–542 (2015).
  * (20) M. Kues, C. Reimer, P. Roztocki, L. R. Cortés, S. Sciara, B. Wetzel, Y. Zhang, A. Cino, S. T. Chu, B. E. Little, D. J. Moss, L. Caspani, J. Azaña, and R. Morandotti, “On-chip generation of high-dimensional entangled quantum states and their coherent control,” Nature 546, 622–626 (2017).
  * (21) R.-B. Jin, R. Shimizu, M. Fujiwara, M. Takeoka, R. Wakabayashi, T. Yamashita, S. Miki, H. Terai, T. Gerrits, and M. Sasaki, “Simple method of generating and distributing frequency-entangled qudits,” Quantum Sci. Technol. 1, 015004 (2016).
  * (22) R.-B. Jin, K. Tazawa, N. Asamura, M. Yabuno, S. Miki, F. China, H. Terai, K. Minoshima, and R. Shimizu, “Quantum optical synthesis in 2D time–frequency space,” APL Photon. 6, 086104 (2021).
  * (23) J. Li, J. Su, L. Cui, T. Xie, Z. Y. Ou, and X. Li, “Generation of pure-state single photons with high heralding efficiency by using a three-stage nonlinear interferometer,” Appl. Phys. Lett. 116, 204002 (2020).
  * (24) C. L. Morrison, F. Graffitti, P. Barrow, A. Pickston, J. Ho, and A. Fedrizzi, “Frequency-bin entanglement from domain-engineered down-conversion,” APL Photon. 7, 066102 (2022).
  * (25) Y. Yang, X. Xiang, F. Hou, R. Quan, B. Li, W. Li, N. Zhu, T. Liu, S. Zhang, R. Dong, and M. Li, “Inherent resolution limit on nonlocal wavelength-to-time mapping with entangled photon pairs,” Opt. Express 28, 7488–7497 (2020).
  * (26) S. Liu, Y. Cui, E. Karimi, and B. A. Malomed, “On-demand harnessing of photonic soliton molecules,” Optica 9, 240–250 (2022).
  * (27) P. J. Mosley, J. S. Lundeen, B. J. Smith, P. Wasylczyk, A. B. U’Ren, C. Silberhorn, and I. A. Walmsley, “Heralded generation of ultrafast single photons in pure quantum states,” Phys. Rev. Lett. 100, 133601 (2008).
  * (28) N. Cai, W.-H. Cai, S. Wang, F. Li, R. Shimizu, and R.-B. Jin, “Broadband-laser-diode pumped periodically poled potassium titanyl phosphate-sagnac polarization-entangled photon source,” J. Opt. Soc. Am. B 39, 77–82 (2022).
  * (29) S. Liu, Y. Zhang, C. Yang, S. Liu, Z. Ge, Y. Li, Y. Li, Z. Zhou, G. Guo, and B. Shi, “Increasing two-photon entangled dimensions by shaping input-beam profiles,” Phys. Rev. A 101, 052324 (2020).
  * (30) B. Brecht, D. V. Reddy, C. Silberhorn, and M. G. Raymer, “Photon temporal modes: A complete framework for quantum information science,” Phys. Rev. X 5, 041017 (2015).

## Spatial-spectral mapping to prepare the frequency entangled qudits:
supplemental document

In this supplemental document, we present the simulations of joint spectral
intensity (JSI) and the spectra of biphotons by adjusting the width and the
center position of the mask. Figure S1 shows the case of changing the width of
a single slit.

Figure S1: Simulation of the spatial-spectral mapping to prepare high-
dimensional entangled states. (a1-a5): The cross section of the pump beam,
spatially filtered by a single slit with different widths. The beam is divided
into 41 bins, and the slit widths in (a1-a5) are 0, 3, 5, 9, 13 bins
respectively. (b1-b5): The calculated JSIs under the corresponding spatial
masks. (c1-c5): The spectra of the signal photon are obtained by projecting
the JSI onto the vertical axis.

Figure S2 shows the case of changing the position of a single slit (moving
up).

Figure S2: Simulation of the spatial-spectral mapping to prepare high-
dimensional entangled states. (a1-a5): The cross section of the pump beam,
spatially filtered by a single slit at different positions. The beam is
divided into 41 bins; the width of the slit is 5 bins; and the center
positions of the slit in (a1-a5) are 21, 22, 23, 24, and 25, respectively.
(b1-b5): The calculated JSIs under the corresponding spatial masks. (c1-c5):
The spectra of the signal photon are obtained by projecting the JSI onto the
vertical axis.

Figure S3 shows the case of changing the position of a single slit (moving
down).

Figure S3: Simulation of the spatial-spectral mapping to prepare high-
dimensional entangled states. (a1-a5): The cross section of the pump beam,
spatially filtered by a single slit at different positions. The beam is
divided into 41 bins; the width of the slit is 5 bins; and the center
positions of the slit in (a1-a5) are 21, 20, 19, 18, and 17, respectively.
(b1-b5): The calculated JSIs under the corresponding spatial masks. (c1-c5):
The spectra of the signal photon are obtained by projecting the JSI onto the
vertical axis.

Figure S4 shows the case of changing the width and position of two slits.

Figure S4: Simulation of the spatial-spectral mapping to prepare high-
dimensional entangled states. (a1-a5): The cross section of the pump beam,
spatially filtered by two slits with different widths and positions. The beam
is divided into 41 bins; the slit-beam-slit widths are 4-1-4, 5-1-5, 6-1-6,
7-1-7, 8-1-8, respectively. (b1-b5): The calculated JSIs under the
corresponding spatial masks. (c1-c5): The spectra of the signal photon are
obtained by projecting the JSI onto the vertical axis.

Figure S5 shows the case of changing the width and position of three or four
slits.

Figure S5: Simulation of the spatial-spectral mapping to prepare high-
dimensional entangled states. (a1-a5): The cross section of the pump beam,
spatially filtered by three or four slits with different widths and positions.
In (a1, a2, a3), the beam is divided into 101 bins; the slit-beam-slit-beam-
slit widths are 11-2-11-2-11, 13-2-13-2-13, 15-2-15-2-15, respectively. In
(a4), the beam is divided into 41 bins; the slit-beam-slit-beam-slit-beam-slit
width is 4-1-4-1-4-1-4. In (a5), the beam is divided into 101 bins; the slit-
beam-slit-beam-slit-beam-slit width is 12-2-12-2-12-2-12. (b1-b5): The
calculated JSIs under the corresponding spatial masks. (c1-c5): The spectra of
the signal photon are obtained by projecting the JSI onto the vertical axis.

Figure S6 shows the case of utilizing a pump beam with a flat-top distribution
rather than a Gaussian distribution.

Figure S6:  Simulation of the spatial-spectral mapping to prepare high-
dimensional entangled states. (a1-a5): The cross section of the pump beam with
a flat-top distribution, spatially filtered by different slits. The beam is
divided into 41 bins; the slit and beam widths are 41, 14-13-14, 9-7-9-7-9,
5-7-5-7-5-7-5, 2-8-2-8-2-8-2-7-2, respectively. (b1-b5): The calculated JSIs
under the corresponding spatial masks. (c1-c5): The spectra of the signal
photon are obtained by projecting the JSI onto the vertical axis.


# Frequency Estimation under Local Differential Privacy [Experiments, Analysis
and Benchmarks]

Graham Cormode, Samuel Maddock, Carsten Maple  University of Warwick, UK
<EMAIL_ADDRESS>
###### Abstract.

Private collection of statistics from a large distributed population is an
important problem, and has led to large scale deployments from several leading
technology companies. The dominant approach requires each user to randomly
perturb their input, leading to guarantees in the local differential privacy
model. In this paper, we place the various approaches that have been suggested
into a common framework, and perform an extensive series of experiments to
understand the tradeoffs between different implementation choices. Our
conclusion is that for the core problems of frequency estimation and heavy
hitter identification, careful choice of algorithms can lead to very effective
solutions that scale to millions of users.

PVLDB Reference Format:
PVLDB, 14(1): XXX-XXX, 2020.
doi:XX.XX/XXX.XX ††This work is licensed under the Creative Commons BY-NC-ND
4.0 International License. Visit https://creativecommons.org/licenses/by-nc-
nd/4.0/ to view a copy of this license. For any use beyond those covered by
this license, obtain permission by emailing<EMAIL_ADDRESS>Copyright is held
by the owner/author(s). Publication rights licensed to the VLDB Endowment.
Proceedings of the VLDB Endowment, Vol. 14, No. 1 ISSN 2150-8097.
doi:XX.XX/XXX.XX

##  1\. Introduction

Collecting frequency statistics underpins a wide range of modern data
analytics tasks, from simple popularity charts to more complex machine
learning tools. If we can centrally collect and aggregate the raw data, then
producing such statistics involves relatively simple data manipulation.
However, in many cases the inputs can be considered sensitive, describing an
individual’s actions, preferences or characteristics. As such, the data
subjects may be reluctant to divulge their true information, even if the
overall statistics would be valuable and socially useful. For example,
consider a mobile app that has access to a user’s location. Understanding the
overall population distribution, where people like to visit at particular
times, and how they move around, has shown to be very valuable for social
studies, urban planning, and tracking disease spread. However, any
individual’s location can be sensitive, allowing inferences to be made on
their health (if they attend a clinic), religious preference (if they attend a
place of worship), sexuality (if they attend a group or club), and more. For
these reasons, users of the app would understandably be unwilling to share
their location information freely with the app.

In recent years, a number of privacy-preserving solutions to this frequency
estimation question have been developed and deployed by large technology
companies, such as Apple, Google and Microsoft. These deployments differ in
the details, but share some common DNA; they are all based on the notion of
Local Differential Privacy, an instantiation of the statistical model of
differential privacy which is applied by each user directly to their own data.
At the heart of these protocols, each user is asked a simple question about
their data by a data analyst. The user then determines randomly whether to
answer the question truthfully, or to pick some false response, according to a
specified probability distribution. This randomization gives each user
“plausible deniability”, and protects their privacy. To fully specify a
protocol, we need to determine what questions are asked, and what probability
distribution is used over the set of possible answers. We also need to show
how the data analyst can combine the answers to build a picture of the overall
frequency distribution, and understand the accuracy of the results.

In this overview and evaluation paper, we seek to draw together the various
approaches (both theoretical and applied) that have been proposed for this
core problem, and place them in a common framework. We show that protocols can
be understood as combining four “layers”: (1) A basic “frequency oracle” that
allows frequencies to be estimated over a moderate sized discrete set of
possibilities; (2) A “sketch” that can reduce the domain size of
possibilities; (3) A “heavy hitters” method that can find frequent items
within a large domain; (4) Post-processing techniques that enforce some
constraints on the results to improve accuracy. Existing protocols can then be
expressed by specifying a particular choice of algorithm for each of the four
layers.

Having described the choices in the four layers, we go on to perform a
sequence of experiments that demonstrate the different behaviors of different
choices of each component, in terms of their accuracy and scalability (speed)
as a function of privacy level, domain size, and other parameter settings.
This leads us to recommend particular combinations for this task, and to
deprecate others. Although not our main objective, we suggest some alternative
choices that have not previously been considered (based on sampling hash
functions), and show that these equal best-in-class performance. Meanwhile, we
show that the first approaches that were introduced for these problems and
popularized them are dominated by more advanced techniques that came
afterwards.

##### Outline.

We proceed as follows: Section 2 introduces the problem definitions and key
technical concepts. We describe the four layers of a private heavy hitters
protocol in turn: first we summarize the “frequency oracles” which are the
basis of all solutions in Section 3. The other three layers of a scheme are
domain size reduction (Section 4), heavy hitter search (Section 5), and post-
processing (Section 6). We structure our experimental study into two parts:
calibration and study of the basic frequency estimation task (Section 7), and
comparative evaluation on the heavy hitters task (Section 8). We summarize our
findings in Section 9.

##  2\. Preliminaries

###  2.1. Local Differential Privacy (LDP)

The definition of differential privacy for a randomized procedure requires
that the output distribution for inputs that are similar (“neighboring”). The
additional requirement to satisfy local differential privacy is that this
should apply to the output distribution of each data owner, rather than just
on the overall output of a centralized procedure. More formally, we consider a
collection of $n$ data owners, who each holds some input $x_{i}$. We assume
that each item $x_{i}\in\mathcal{D}$ is drawn from some domain $\mathcal{D}$
of size $|\mathcal{D}|=d$. We write $e^{(j)}$ as the basis vector of dimension
$d$ such that $e^{(j)}_{j}=1$, and 0 otherwise. Each user will apply a local
randomizer $\mathcal{R}$ to their input to obtain an output from a domain
$\mathcal{Y}$, so that $\mathcal{R}:\mathcal{D}\rightarrow\mathcal{Y}$.

###### Definition 2.1 (Local differential privacy).

The local randomized $\mathcal{R}$ is $\epsilon$-locally differentially
private ($\epsilon$-LDP) if for all pairs of inputs
$x\in\mathcal{D},x^{\prime}\in\mathcal{D}$ and outputs $y\in\mathcal{Y}$ we
have that $\frac{\Pr[\mathcal{R}(x)=y]}{\Pr[\mathcal{R}(x^{\prime})=y]}\leq
e^{\epsilon}$.

The objective is then to obtain the best accuracy for a particular task while
guaranteeing the $\epsilon$-LDP property. Further background and discussion of
local differential privacy is presented in recent surveys and tutorials (Yang
et al., 2020; Wang et al., 2020; Cormode et al., 2018).

###  2.2. Frequency Estimation and Heavy Hitters

With the same setup as the previous section, we define the (raw) frequency of
item $x\in\mathcal{D}$ to be $f(x)=|\\{i:x_{i}=x\\}|$. That is, we simply
count the number of occurrences of item $x$ among the population. Natural
variations consider a weighted version of this definition, but our focus here
is on the unweighted (unit weight) notion. Our objective will be to obtain a
(private, randomized) estimator $\hat{f}$ for $f$. We will typically be
interested in estimators $\hat{f}$ that are _unbiased_ , that is
$\mathsf{E}[\hat{f}(x)]=f(x)$, and which have small _error_
$|\hat{f}(x)-f(x)|$. When the estimator is unbiased, the mean squared error is
equivalent to the variance, $\mathsf{E}[(\hat{f}(x)-f(x))^{2}]$. The heavy
hitters are items $x$ whose frequency is large – this is often formalized to
be those items which exceed a threshold of the total weight or total sum of
squares, or are in the top-$k$ largest weights.

###  2.3. Hadamard Transformation

Several private mechanisms depend on the Hadamard transformation, which is a
particular instance of a discrete Fourier transform. It is an orthogonal
transform $\phi$, described by a $d\times d$ matrix with
$\phi_{i,j}=d^{-1/2}(-1)^{\langle i,j\rangle}$, where $d$ is a power of 2, and
$\langle i,j\rangle$ counts the number of indices in the binary
representations of $i$ and $j$ where both have 1 bit. Given a vector $x$ of
dimension $d$, the Hadamard transformation of $x$ is the vector of
coefficients $\Theta=\phi x^{T}$. For general $x$ this can be computed
efficiently in time $O(d\log d)$, but we often apply this to sparse inputs and
find one coefficient, which can be done with a constant number of mathematical
operations.

##  3\. Frequency Oracles

The notion of a frequency oracle has emerged from a sequence of works in
differential privacy as an LDP mechanism for frequency estimation. Oracles
vary in their construction, accuracy guarantees, and the size of domain for
which they are best suited.

###  3.1. “Pure” Protocols

A useful abstraction introduced by Wang et al. (Wang et al., 2017) is the
notion of a ‘pure’ protocol for frequency estimation. A protocol is considered
to be pure if its output can be interpreted as expressing “support” for a
subset of possible input values, encoded as $\operatorname{supp}(y)$. Pure
protocols are characterized by two probabilities, $p^{*}$ and $q^{*}$, so that
$\Pr[\mathcal{R}(x)\in\\{y:x\in\operatorname{supp}(y)\\}]=p^{*}$, which
encodes the probability that input $x$ is mapped to an output that supports
$x$; and
$\Pr[\mathcal{R}(x^{\prime})\in\\{y:x\in\operatorname{supp}(y)\\}]=q^{*}\text{
for all $x^{\prime}\neq x$}$, which encodes the (lower) probability that an
input $x^{\prime}$ is mapped to an input that supports $x$ for every
$x^{\prime}\neq x$.

Many of the frequency oracles below are instances of a pure protocol, which
makes their analysis convenient. The $\epsilon$-LDP condition is satisfied
provided $p^{*}/q^{*}\leq e^{\epsilon}$. To aggregate the responses from a
pure protocol to make an estimate $\hat{f}(x)$, for each report
$y_{i}=\mathcal{R}(x_{i})$, we compute
$\hat{f}_{i}(x)=\frac{\mathcal{I}_{x\in\operatorname{supp}(y_{i})}-q^{*}}{p^{*}-q^{*}}$
where $\mathcal{I}_{P}$ is the function that is 1 if $P$ is true, and 0
otherwise. We can verify that

$\mathsf{E}[\hat{f}_{i}(x)]={(\mathsf{E}[\mathcal{R}(x_{i})\in\\{y:x\in\operatorname{supp}(y)\\}]-q^{*})}/{(p^{*}-q^{*})}$

which simplifies to 1 if $x_{i}=x$, and 0 if $x_{i}\neq x$.

The contribution to the (conditional) variance of this estimator for each user
$i$ can be computed for the two cases, as:

(1)
$\displaystyle\mathsf{Var}[\hat{f}_{i}(x)]=\frac{\mathsf{Var}[\mathcal{I}_{x\in\operatorname{supp}(y_{i})}]}{(p^{*}-q^{*})^{2}}=$
$\displaystyle\frac{p^{*}(1-p^{*})}{(p^{*}-q^{*})^{2}}\\!\\!\\!\\!\\!\\!\\!\\!$
if$x=x_{i}$ (2)
$\displaystyle\mathsf{Var}[\hat{f}_{i}(x)]=\frac{\mathsf{Var}[\mathcal{I}_{x\in\operatorname{supp}(y_{i})}]}{(p^{*}-q^{*})^{2}}=$
$\displaystyle\frac{q^{*}(1-q^{*})}{(p^{*}-q^{*})^{2}}\\!\\!\\!\\!\\!\\!\\!\\!$
if$x\neq\\!x_{i}$

The overall estimate $\hat{f}(x)$ is the sum of the individual estimates, and
its variance is the sum of these variances. The exact value of the overall
estimate’s variance will depend on the true frequency of $x$ in the input. If
we assume that the domain is large and no input dominates the others, then
this expression is determined by the expression in $q$, and we can focus on
the contribution to the variance from the $x\neq x_{i}$ case in (2).

###  3.2. Direct Encoding (DE)

The idea of direct encoding (also known as randomized response) is to set the
output range to be the same as the input, i.e., $\mathcal{Y}=\mathcal{D}$, and
to have some probability of providing a “false” response, while maximizing the
likelihood of giving a “truthful” response (Erlingsson et al., 2014; Wang et
al., 2017). If we set the truthful probability to be
$\Pr[\mathcal{R}(x)=x]=p$, then for any $y\neq x$ we must have
$\Pr[\mathcal{R}(x)=y]=(1-p)/(d-1)=q$ and LDP requires $p/((1-p)/(d-1))\leq
e^{\epsilon}$. Equating and rearranging this, we obtain
$p=\frac{e^{\epsilon}}{e^{\epsilon}+d-1}$ and $q=\frac{1}{e^{\epsilon}+d-1}$.
This yields a pure protocol with $p^{*}=p$ and $q^{*}=q$. To build an unbiased
estimator for this mechanism, the contribution from each report is
$\hat{f}_{i}(x)=\frac{(e^{\epsilon}+d-1)\mathbb{I}_{x}(\mathcal{R}(x_{i}))-1}{e^{\epsilon}-1}$,
where $\mathbb{I}_{x}(\cdot)$ is the indicator function for whether its
argument is equal to $x$. It then follows that
$\mathsf{E}[\hat{f}_{i}(x)]=\mathbb{I}_{x}(\mathcal{R}(x_{i}))$, and that

$\mathsf{Var}[\hat{f}_{i}(x)]=(\mathbb{I}_{x}(\mathcal{R}(x_{i}))p(1-p)+(1-\mathbb{I}_{x}(\mathcal{R}(x_{i}))q(1-q))\left(\frac{e^{\epsilon}+d-1}{e^{\epsilon}-1}\right)^{2}.$

The variance term corresponding to case $x_{i}\neq x$ from (2) for DE is
$\mathsf{Var}[\hat{f}(x)]=\frac{q(1-q)}{(p-q)^{2}}=\frac{e^{\epsilon}+d-2}{(e^{\epsilon}-1)^{2}}.$

###  3.3. Unary Encodings (OUE, SUE)

Unary encoding methods interpret the user’s input $x_{i}$ as a one-hot
$d$-dimensional vector $e^{(x_{i})}$, and independently perturb each bit.
Hence, $\mathcal{Y}=\\{0,1\\}^{d}$. Suppose we retain the (sole) input 1 with
probability $p$, and flip each 0 to a 1 with probability $q$. This provides a
pure protocol with $p=p^{*}$ and $q=q^{*}$, where we have the freedom to
choose $p$ and $q$ subject to satisfying $p(1-q)/(1-p)q\leq e^{\epsilon}$ to
ensure $\epsilon$-LDP111Considering neighboring inputs which differ in the
location of the solitary 1 bit.. “Symmetric” unary encoding (SUE), as used in
RAPPOR (Erlingsson et al., 2014) sets $p+q=1$, yielding
$p=\frac{e^{\epsilon/2}}{e^{\epsilon/2}+1}$ and
$q=\frac{1}{e^{\epsilon/2}+1}$. The corresponding variance is
$pq/(p-q)^{2}=e^{\epsilon/2}/(e^{\epsilon/2}-1)^{2}$. “Optimized” unary
encoding (OUE, (Wang et al., 2017)) aims to minimize (2) subject to
$p(1-q)/(1-p)q=e^{\epsilon}$. Rearranging this sets
$p=q\frac{e^{\epsilon}}{1+qe^{\epsilon}-q}$, and hence
$p-q=\frac{(e^{\epsilon}-1)q(1-q)}{1+qe^{\epsilon}-q}$ and so (2) gives the
variance contribution as
$\frac{(1+q(e^{\epsilon}-1))^{2}}{(e^{\epsilon}-1)^{2}q(1-q)}$. We can
minimize the factor in $q$ by taking the partial derivative with respect to
$q$ and equating to zero. This leads us to choose $q=1/(e^{\epsilon}+1)$, and
the expression for (2) simplifies to
$\frac{4e^{\epsilon}}{(e^{\epsilon}-1)^{2}}$. This fixes $p=\frac{1}{2}$,
i.e., the encoding retains the sole original input 1 with constant
probability, while flipping zeros only with smaller probability
$1/(e^{\epsilon}+1)$.

###  3.4. Local Hashing (BLH, OLH, FLH)

As $d$ gets larger, the dependence on $d$ for DE is undesirable, and
subsequent approaches remove this dependency. The local hashing approach has
each user pick a hash function with which to encode their input, and send back
the hash function along with the hashed input encoded via direct encoding
(Wang et al., 2017). Here, the domain of the hash is chosen to be $g$, with
$g$ chosen as a function of $\epsilon$, typically much smaller than $d$. That
is, user $i$ picks a (universal) hash function $h_{i}$ that maps
$[d]\rightarrow[g]$. Next, the user applies direct encoding to the hashed
value $h_{i}(x_{i})$, and sends this output along with the description of
$h_{i}$. From this, the aggregator can obtain an unbiased estimator for
$\hat{f}(h_{i}(x))$. This gives a pure protocol, where the probability that
$x_{i}$ is mapped to the hash value that supports $x_{i}$ is
$p^{*}=\frac{e^{\epsilon}}{e^{\epsilon}+g-1}$, while the probability that
$x\neq x_{i}$ is mapped to that hash value becomes the uniform chance
$q^{*}=1/g$.

The variance contribution due to (2)
$\frac{(e^{\epsilon}+g-1)^{2}}{(g-1)(e^{\epsilon}-1)^{2}}$. A simple case is
Binary Local Hashing (BLH), which just sets $g=2$, and obtains variance
$\frac{(e^{\epsilon}+1)^{2}}{(e^{\epsilon}-1)^{2}}$. However, by choosing
$g-1=e^{\epsilon}$ (Optimal Local Hashing, OLH) we minimize the factor
$\frac{(e^{\epsilon}+(g-1))^{2}}{(g-1)}$, and hence obtain a variance of
$\frac{4e^{\epsilon}}{(e^{\epsilon}-1)^{2}}$. Note that this leads us to
choose $p^{*}=\frac{1}{2}$, i.e., we report the true hash value with
probability half. A limitation of the OLH approach is that it can be very
slow. This is since we need to make $O(nd)$ hash function calls to calculate
which domain elements a user’s perturbed item contributes frequency towards.

Here, we propose a “fast” heuristic modification which aims to speed up the
procedure, dubbed Fast Local Hashing (FLH). On the client-side, instead of
sampling a hash function uniformly at random from some universal hash family,
we introduce a new parameter $k^{\prime}$ and restrict clients to uniformly
choosing from $k^{\prime}$ hash functions. Hence, we sacrifice some
theoretical guarantees on accuracy in order to achieve computational gains on
the server-side aggregation. In practice, it turns out that even for small
values of $k^{\prime}$, FLH performs reasonably well.

The protocol is largely unchanged from OLH. On the client-side, we sample a
hash function uniformly at random from $\\{h_{1},\dots h_{k^{\prime}}\\}$. On
the server-side, we pre-compute a $k^{\prime}\times d$ matrix where each row
corresponds to a hash function and each column refers to a domain value. We
set the $(i,j)^{th}$ entry to $h_{i}(j)$ to reduce the total number of hash
function calls from $O(nd)$ to $O(k^{\prime}d)$. The server collates all
reports for the same hash function, and process them in a batch. To gain a
speedup over OLH, we must choose $k^{\prime}\ll n$. In our experiments
(Section 7.1), we empirically validate FLH, and show that it approaches OLH in
accuracy when only sampling a small number of hash functions, achieving a
substantial speed up.

###  3.5. Hadamard Encodings (HM, HR)

We next consider two related approaches based on the Hadamard transform. The
Hadamard mechanism (HM) samples one (or more) coefficients from the Hadamard
transform of the users input (treated as a sparse vector), and applies direct
encoding to the result (app, 2017; Bassily et al., 2020). So
$\mathcal{Y}=\\{-1,+1\\}$ and given a user’s input $x_{i}\in\mathcal{D}=[d]$,
we sample an index $j\in[d]$, and compute the (scaled-up) coefficient
$\theta_{j}^{(i)}=\phi_{x_{i},j}=(-1)^{\langle x_{i},j\rangle}$ (see Section
2.3). Then with probability $p=\frac{e^{\epsilon}}{1+e^{\epsilon}}$, the
mechanism reports $(j,\theta_{j}^{(i)})$, otherwise it reports
$(j,-\theta_{j}^{(i)})$. Since
$\frac{p}{1-p}=\frac{e^{\epsilon}}{1+e^{\epsilon}}\cdot\frac{1+e^{\epsilon}}{1}=e^{\epsilon}$,
this meets the LDP guarantee.

An advantage of HM is that it is much faster to process reports from users,
since the (fast) Hadamard inverse transform can be applied to the (aggregated)
reports. To build an unbiased estimator for frequencies, we take the
contribution of the inverse of the unbiased estimator of the reported
$\hat{\theta}_{j}=\sum_{i}\theta_{j}^{(i)}$. The unbiased estimator for
$\theta_{j}$ is ${\hat{\theta}_{j}}/{(2p-1)}$. For a given $x$, to estimate
$f(x)$, we sum the contribution to $f(x)$ from all reports. That is,
$\hat{f}_{i}(x)=\sum_{j}\phi_{x,j}\cdot{\hat{\theta}_{j}^{(i)}}/{(2p-1)}$ and
$\hat{f}(x)=\sum_{i}\hat{f}_{i}(x)$ as before. We can interpret the $j$’th row
of the Hadamard matrix as defining a (universal) hash function on the input
domain with range $g=2$. The analysis for BLH applies, and the variance bound
(2) is $\frac{(e^{\epsilon}+1)^{2}}{(e^{\epsilon}-1)^{2}}$.

To improve this bound when $e^{\epsilon}$ is large, we can sample $t$ Hadamard
coefficients, to produce a hash function with $g=2^{t}$ possible outcomes.
This preserves the result with probability
$p^{*}=\frac{e^{\epsilon}}{e^{\epsilon}+2^{t}-1}$, and otherwise perturbs it
uniformly. The resulting variance from (2) is
$\frac{(e^{\epsilon}+2^{t}-1)^{2}}{(2^{t}-1)(e^{\epsilon}-1)^{2}}$. If we can
choose $t$ so that $2^{t}-1=e^{\epsilon}$,222This is possible when
$e^{\epsilon}$ is moderately large, i.e., $e^{\epsilon}=3,7,15.\ldots$ then we
obtain the optimal variance of $\frac{4e^{\epsilon}}{(e^{\epsilon}-1)^{2}}$.
We can aggregate quickly by interpreting the hash value as $t$ Hadamard
coefficients, each of which is transmitted correctly with probability
$\frac{e^{\epsilon}}{e^{\epsilon}+2^{t}-1}+\frac{2^{t-1}-1}{e^{\epsilon}+2^{t}-1}\cdot\frac{2^{t-1}-1}{2^{t}-1}=\frac{e^{\epsilon}+2^{t-1}-1}{e^{\epsilon}+2^{t}-1}$.
For $e^{\epsilon}=3$ and $t=2$, this fraction is $2/3$.

The Hadamard response protocol, HR, (Acharya et al., 2019) also makes use of
the Hadamard matrix, and is very similar to HM. Each user with input $x_{i}$
will report the index of some Hadamard coefficient, i.e., some $j$, where the
corresponding coefficient value is $\theta_{j}=\phi(e^{(x_{i})})^{T}$.
Equivalently, the set of all Hadamard coefficient values $\Theta$ correspond
to the $x_{i}$’th column of Hadamard matrix $\phi$. With probability
$p=\frac{e^{\epsilon}}{e^{\epsilon}+1}$, we report a $j$ such that
$\theta_{j}=+1$, otherwise (with probability $q=\frac{1}{e^{\epsilon}+1}$), we
report a coefficient $j$ with $\theta_{j}=-1$. This again represents a pure
protocol, since by properties of the Hadamard transform, $p^{*}=p$, while
$q^{*}=\frac{1}{2}$. Hence, we compute our variance bound as
$(\frac{e^{\epsilon}+1}{e^{\epsilon}-1})^{2}$. For higher privacy (smaller
epsilon), a variant protocol is based on the Hadamard matrix, but is no longer
pure, since the probability of $x^{\prime}$ being mapped to the support of $x$
is no longer the same for all $x^{\prime}$. This construction achieves a
variance bound of $O(\frac{e^{\epsilon}}{(e^{\epsilon}-1)^{2}})$.

###  3.6. Summary of Frequency Oracles

Table 1. Comparison of Frequency Oracle costs FO | Variance bound | Decode time | Communication
---|---|---|---
DE | $O((e^{\epsilon}+d-2)/(e^{\epsilon}-1)^{2})$ | $O(n+d)$ | $O(\log d)$
SUE | $O(e^{\epsilon/2}/(e^{\epsilon/2}-1)^{2})$ | $O(nd)$ | $O(d)$
OUE | $O(e^{\epsilon}/(e^{\epsilon}-1)^{2})$ | $O(nd)$ | $O(d)$
BLH | $O((e^{\epsilon}+1)^{2}/(e^{\epsilon}-1)^{2})$ | $O(nd)$ | $O(\log d)$
OLH | $O(e^{\epsilon}/(e^{\epsilon}-1)^{2})$ | $O(nd)$ | $O(\log d)$
FLH | (heuristic) | $O(k^{\prime}d)$ | $O(\log d)$
HM | $O((e^{\epsilon}+1)^{2}/(e^{\epsilon}-1)^{2})$ | $O(n+d)$ | $O(\log d)$
HR | $O(e^{\epsilon}/(e^{\epsilon}-1)^{2})$ | $O(n+d)$ | $O(\log d)$

The concept of a frequency oracle is an abstraction of a private protocol
allowing us to accurately estimate the frequency of items from a domain. The
listed approaches achieve different combinations of accuracy, speed and
communication costs in order to guarantee $\epsilon$-LDP. We summarize the key
properties in Table 1.

Across all frequency oracles, the variance obtained is quite similar, with a
dependence on $\frac{1}{(e^{\epsilon}-1)^{2}}$. When $\epsilon$ is suitably
small, we have $e^{\epsilon}\approx 1+\epsilon$, and so
${(e^{\epsilon}-1)^{-2}}\approx{\epsilon^{-2}}$. The numerator of the variance
term is typically either $(e^{\epsilon}+1)^{2}$ (for basic Hadamard
mechanism), or $4e^{\epsilon}$ (for Hadamard response, optimized local hashing
and optimized unary encoding). Since $(x+1)^{2}-4x=(x-1)^{2}>0$ for $x>1$, so
$4e^{\epsilon}<(e^{\epsilon}+1)^{2}$, and hence we would expect these last
methods to obtain better accuracy, as we test in Section 7.

##  4\. Domain size reduction

When working with data over large domain sizes $d$, performing a search for
frequent items can be slow, and incur more false positives as we query more
possibilities. A natural approach is to encode the input in a smaller domain,
via techniques such as “sketching” (Cormode and Yi, 2020). We can first encode
each user’s input into a small sketch, then apply the frequency oracle to each
sketch, so that an aggregator can reconstruct a sketch to query. Since the
oracles described above all rely on the user’s input being sparse (i.e.,
encoded by a 1-hot vector), we seek sketches which are sparsity-preserving.
Each sketch is defined by an array with $r$ rows and $c$ columns, and some
hash functions which map items to the array. Next, we survey three sketch
approaches that have been used for this purpose.

Bloom Filter. The Bloom filter (BF) is a well-known way to represent a set of
items (Bloom, 1970). It makes use of $k$ hash functions $h_{1}\ldots h_{k}$
which map $[d]\rightarrow[m]$. This generates a one-dimensional binary array
of length $m$, whose $j$’th entry is 1 if there is a hash function
$h_{\ell}(x)=j$. Encoding a single user’s item $x_{i}$ means that we obtain an
encoding with at most $k$ 1’s set. Frequency oracles can be applied by first
sampling one of the hash functions $h_{\ell}$, and then processing the
resulting one-hot encoding of $h_{\ell}(x_{i})$. Decoding the result and
obtaining an estimate from the noisily reconstructed Bloom filter can be done
e.g., using regularized regression as a heuristic (Erlingsson et al., 2014).

Count-Min Sketches. The Count-Min (CM) sketch is similar to the Bloom Filter,
but now hash function $h_{\ell}$ maps to row $\ell$ of a $r\times c$ sized
array (Cormode and Muthukrishnan, 2005). Item $x$ is then mapped to entry
$[\ell,h_{\ell}(x)]$ for each $\ell$, so the sketch can be treated as $r$ one-
hot vectors, one for each of the $r$ hash functions. The user first samples
$\ell$ uniformly from $\\{1\ldots r\\}$, and encodes $h_{\ell}(x_{i})$ using a
frequency oracle. The aggregator reconstructs a noisy sketch, and can estimate
$f(x)$ by considering the weights associated with $h_{\ell}(x)$: this can be
via taking the minimum, the mean, or the median (giving the Count-Mean or
Count-Median sketch, respectively, with appropriate corrections to make the
estimation unbiased (app, 2017; Cormode and Yi, 2020)). The variance
introduced due to sketching depends on $c$, the length of each row. Smaller
$c$ and $r$ reduces the cost of storing and analyzing the sketch, but larger
$c$ improves the accuracy.

Count Sketches. The Count sketch (CS) follows the description of the CM
sketch, but with one further twist: another hash function $g_{\ell}$
determines whether to represent $x$ with a $+1$ or $-1$ value (Charikar et
al., 2002). That is, the user samples $\ell$ uniformly from $r$, and encodes
$h_{\ell}(x_{i})$ with a weight of $g_{\ell}(x_{i})$, where
$g_{\ell}(x_{i})=\\{-1,+1\\}$. This is most easily handled by the Hadamard
mechanism, which handles negative inputs without any modification. Other FOs
can be extended to handle negative weights, such as by handling positive and
negative weights separately. To estimate $f(x)$, the aggregator takes the
weight associated with $h_{\ell}(x)$, multiplied by $g_{\ell}(x)$. The
variance due to sketching also depends on $1/c$, as in the previous case.

##  5\. Heavy Hitter Search

Domain size reduction reduces the size of the objects manipulated in frequency
estimation protocols. However, if the original domain within which the items
reside is very large, it is still costly for the data analyst to search
through the set of all possibilities to find those with high frequency. For
example, the input items may be long strings (such as URLs), and it would be
too slow to enumerate all strings up to a fixed length. This is known as the
“heavy hitters” problem, and has been addressed extensively in the (non-
private) data streaming setting. Approaches that work in the local privacy
model are based on ideas from coding theory and group testing.

Sequence Fragment Puzzle (SFP). The core idea of the “sequence fragment
puzzle” (SFP) approach is that for a long string that is frequent within a
population, every substring will also be (at least as) frequent (app, 2017).
Naively, we might try to find all frequent substrings and glue them back
together. However, this risks creating false positives. For example, suppose
“apples” and “orange” are both frequent (6 character) strings. We could break
these into 3 character substrings, as “app”, “les”, “ora” and “nge”, which
will all be frequent. But when reconstructing, we should avoid reporting non-
frequent strings, such as “orales” or “ngeapp”. A heuristic approach is for
each user to provide information on pairs of disjoint substrings that co-occur
together in their string – e.g., dividing strings into four substrings so that
one user may reporting on their first and third substring, another on their
third and fourth. This gives the aggregator partial information about which
substrings co-occur frequently, allowing statistical inference to rebuild the
whole strings (Fanti et al., 2016).

The SFP method gives a principled solution with two extra steps. The first is
to provide the index of each substring within the string, and the second is to
also tag each substring with a hash of the whole string. In our example, the
pieces might be represented as (1, “app”, hash(apples)), (2, “les”,
hash(apples)), (1, “ora”, hash(orange)), (2, “nge”, hash(orange)). If an
analyst detects these four tuples as frequent then, provided there are no hash
collisions, the only way to combine them into strings yields the desired
“apples” and “orange”. In the LDP setting, each user takes their string and
chooses a substring, then encodes this substring concatenated with the hash
value via a frequency oracle (possibly with domain size reduction). The index
of the substring not need be kept private, and can be sent without
perturbation. This allows the analyst to find which substrings are globally
frequent by enumerating all possibilities, and put these together to solve the
‘puzzle’.

There are several parameters to choose for SFP: (1) The length of the
substrings to consider. Too short, and there is a lot of work to recombine the
strings; too long, and the time cost of the analyst to search the frequency
oracles is too high. Apple’s implementation uses substrings of length 2-3
characters over the Roman alphabet. (2) The size of the hash function. This
should be large enough to avoid collisions between frequent strings, but not
so large that it makes the concatenated message large. Apple’s implementation
uses an 8-bit hash. (3) How each user should select a substring to report. The
suggested approach is to pad each user string to the same length, and
uniformly sample a single index to report. In Apple’s example, strings are 10
characters long, and we sample one of the five 2-character strings starting at
an odd location.

Reporting a single substring per user achieves better accuracy than reporting
multiple substrings with a shared privacy budget. To avoid false positives, we
can additionally maintain a frequency oracle over the full strings (at some
privacy cost).

Hierarchical Search (PEM, TH). An alternative approach to finding heavy
hitters is hierarchically based. Considering the inputs as strings, observe
that if a string is a heavy hitter, then every prefix will be at least as
frequent within the population. So we can build frequency oracles on prefixes
of increasing length, and use these to step towards the full length strings.
The key to making this efficient is the converse of the above statement: if a
prefix is not a heavy hitter, then no extension of the prefix can be a heavy
hitter either.

The idea of hierarchical methods is to first enumerate all prefixes of some
fixed length (say, 2 or 3 alphabet characters), and estimate their frequency
with an appropriate frequency oracle. This yields a subset of heavy hitter
prefixes. The next phase uses a frequency oracle on prefixes of a longer
length, and estimates the frequency of all possible extensions of the heavy
hitter prefixes. This is repeated until we reach the limiting length of
prefixes for which a frequency oracle has been built. Variants on this idea
are known under various names – prefix extending method (PEM) (Wang et al.,
2019), TreeHistogram (TH) (Bassily et al., 2020), and PrivTrie (Wang et al.,
2018), but we group them together under the hierarchical search banner. The
chief difference is how they report information: in PEM, each user is equally
likely to report each prefix length, while in TH the privacy budget is split
evenly to report on both the full string and and a uniformly chosen shorter
prefix. As with SFP, there are several parameters to choose: (1) The initial
prefix length, and how long to extend the prefix by in each phase; and (2) How
each user should decide what information to report. A natural choice is for
each user to uniformly select a prefix length from those specified, and only
encode their input once.

Error-correcting codes. We briefly mention that other approaches have been
proposed based on applying error-correcting codes (ECC) to allow recovering
the heavy hitter strings. This general approach was first presented by Bassily
and Smith (Bassily and Smith, 2015), and used in the Bitstogram protocol
(Bassily et al., 2020) to give improved theoretical bounds. However, we omit
this approach from further consideration for brevity, since it has not been
widely implemented.

##  6\. Post-processing

We finally address how additional post-processing can be done on the outputs
of frequency estimation and heavy hitter protocols with the aim of improving
accuracy. Post-processing is useful, since the output of frequency oracles can
be quite noisy: we can obtain negative frequencies, or outputs which sum to
more the number of inputs. We consider the following options in our
experiments:
— No Normalization: No modification of the estimated frequencies.
— Non-Negative Normalization: We round all negative frequency estimates to
$0$.
— Additive Normalization: We round negative estimates to $0$. For the rest of
the values, we add/subtract some constant $\delta$ to ensure that
$\sum_{v\in\mathcal{D}_{>0}}(\tilde{f}(v)+\delta)=n$ where
$\mathcal{D}_{>0}=\\{v:\tilde{f}(v)>0\\}$.
— Projection onto the probability simplex: Here we project our frequency
oracle onto the probability simplex, ensuring that
$\sum_{v\in\mathcal{D}}\tilde{f}(v)=n,\tilde{f}(v)\geq 0$. See (Acharya et
al., 2018) and (Wang and Carreira-Perpiñán, 2013) for more information.
— Threshold Cut: For threshold-cut, when estimating the whole domain we sort
our frequency estimates in decreasing order and keep them until we get a total
frequency which is $>n$. At this point we round every remaining estimate down
to $0$.

Note that applying some of these post-processing steps can be quite time
consuming when the domain is large. Non-negative normalization is
straightforward, but other post-processing requires enumerating all possible
items in the domain to determine the normalization parameters. Thus, these
can’t be applied for frequency oracles over very large input domains.

##  7\. Frequency Estimation Experiments

Our first experiments focus on the ability of different frequency oracles and
their variants to estimate frequencies. We are mainly interested in measuring
the mean squared error (MSE) of the frequencies estimated over the whole
domain i.e., $\frac{1}{d}\sum_{x\in\mathcal{D}}(\hat{f}(x)-f(x))^{2}$, which
also serves as an empirical estimation of the variance.

All experiments were implemented in Python 3.7.4 and run on a Windows PC, with
an AMD Ryzen 5 3600 3.6GHz CPU and 16GB of RAM. The code is publicly available
333https://github.com/Samuel-Maddock/pure-LDP. Unless stated otherwise all,
experiments were run 5 times with their results averaged. We test the
frequency oracles mainly on synthetic data generated from a Zipf distribution
with skewness $s=1.1$.

(a)

(b)

(c)

Figure 1. DE, SUE, OUE Experiments

###  7.1. Direct and Unary Encoding Methods

We start by comparing Direct Encoding (DE), Symmetric Unary Encoding (SUE) and
Optimal Unary Encoding (OUE) (Sections 3.2 and 3.3). In Figure
LABEL:fig:group1_eps we see the effects of varying the privacy budget
$\epsilon\in\\{0.5,1,\dots,5\\}$ while fixing $d=1024,n=100{,}000$. Although
the domain is not very large, OUE/SUE both outperform DE due to the high
amount of noise needed in DE. For relatively low values of $\epsilon$ (high
privacy) both OUE and SUE perform quite similarly, however for larger values
of $\epsilon$ we can see OUE has the lower MSE matching the theoretical
results.

In Figure LABEL:fig:group1_d we vary the domain size
$d\in\\{2^{2},2^{3},\dots,2^{11}\\}$ while fixing $\epsilon=3,n=100{,}000$.
For very small domains, DE has much lower variance than OUE/SUE, however at
around $d=2^{6}$ OUE/SUE start to outperform DE. This agrees with the
theoretical results, since the variance of DE scales in $d$ as we require
significantly more noise in order to maintain LDP across the domain.

In practice DE is fast and simple to implement, since all that needs to be
returned by a client is a single value in the domain. OUE/SUE require
perturbing a $d$-length one-hot vector and sending it to the server. When the
domain is large, this is impractical in terms of communication cost. Our
results corroborate those of (Wang et al., 2017).

Figure 2. Effect of varying $k^{\prime}$ in FLH,
$d=500,n=1,\\!000,\\!000,\epsilon=3$

(a)

(b)

(c)

Figure 3. Local Hashing (BLH, OLH, FLH) Experiments

###  7.2. Local Hashing Methods

We compare Optimised Local Hashing (OLH) with Binary Local Hashing (BLH)
(i.e., using a hash function with $g=2$ possible values) and with Fast Local
Hashing (FLH). The parameter $k^{\prime}$ in FLH determines the number of hash
functions used, with larger $k^{\prime}$ resulting in more accurate estimates
at the expense of slower server-side aggregation. In Figure 2 we fix
$d=500,n=1{,}000{,}000$ and vary $k^{\prime}$ between $100$ and $20{,}000$ in
increments of $100$. The blue line represents the MSE of FLH as we vary
$k^{\prime}$, and the red line represents the average OLH error on the same
data. As $k^{\prime}$ increases, the MSE of FLH approaches that of OLH. At
around $k^{\prime}=10{,}000$ the accuracy of FLH is indistinguishable from
that of OLH. When $k^{\prime}=10{,}000$, FLH’s aggregation is $\approx 6$x
faster than OLH.

Figure 3 shows experiments for varying both $\epsilon$ and $d$ for OLH, BLH
and FLH with $k^{\prime}=100,500,1000,5000,10{,}000$. Figure
LABEL:fig:group2_eps shows the effect of varying
$\epsilon\in\\{0.5,1,\dots,5\\}$ with $d=1024,n=100,000$. We observe that the
MSE of BLH is consistent for different $\epsilon$ values, and as noted by
(Wang et al., 2017) this is likely due to the use of a binary hash function
losing a lot of information, making BLH relatively robust to low $\epsilon$
values. OLH has the lowest MSE followed closely by the FLH oracles in
decreasing order of their $k^{\prime}$ value, as we would expect. We see
similar relative behavior as we vary $d$ in Figure LABEL:fig:group2_d.

Meanwhile, in Figure LABEL:fig:group2_time, we plot the total time taken
across all client/server calculations as we vary $d$. Since we are enumerating
the entire domain, this naturally has a linear dependency on $d$.
Nevertheless, we can see the effect that small values of $k$ in FLH have in
terms of reducing the total time taken compared to OLH and BLH, which both
have to iterate over $n$ different hash functions. A smaller $k^{\prime}$
parameter reduces the aggregation time of the server, but both client-side
perturbation and server-side estimation take the same amount of time.

(a)

(b)

(c)

Figure 4. Hadamard Experiments

###  7.3. Hadamard Methods

We compare the two Hadamard-based approaches (Section 3.5): Hadamard Mechanism
(HM) and Hadamard Response (HR). The main parameter of the Hadamard Mechanism
is the number of Hadamard coefficients to sample, which we denote as $t$. In
Figure LABEL:fig:group3_t we fix $d=1024,n=100{,}000$ while varying $t$ along
the $x$-axis with each line representing a different value of $\epsilon$. We
see that for low values of $\epsilon$ (high privacy) that the optimal value is
$t=1$, whereas for $\epsilon=3$, $t=3$ has the lowest MSE. Generally we see
optimal values for $t=\lceil\epsilon\rceil$, consistent with our analysis of
this mechanism.

In Figure LABEL:fig:group3_eps we vary $\epsilon$ while fixing
$d=1024,n=100{,}000$. We compare the MSE of HR with the HM with $t=1,2,3,4,5$.
We can see that for the smaller $\epsilon$ values HR and HM perform very
similarly when using the optimal values of $t=1,2$. Additional experiments
were performed in the low privacy regime ($\epsilon<1$) and we found that the
performance of HR and HM were almost identical but we omit these plots for
brevity. For higher $\epsilon$ values ($\epsilon>3$), HR starts to outperform
the HM oracles, which is quite notable for methods which appear so similar on
paper. In Figure LABEL:fig:group3_d we vary $d$ while fixing
$\epsilon=3,n=100{,}000$. We can again see that HM $t=3$ has the best MSE for
$\epsilon=3$, but still more than that of HR for the large domains. For very
small domains, the behavior is quite variable, but nevertheless HM appears
preferable in this regime.

(a)

(b)

(c)

Figure 5. Comparison Experiments

###  7.4. Cross-comparison of Frequency Oracles

In Figure 5 we compare OUE, OLH, FLH, HR and HM, varying $\epsilon,d$ and
looking at the total time taken. We omit DE since we have seen that its
performance decreases in larger domains, and BLH/SUE since OLH/OUE have better
performance. For FLH we take $k^{\prime}=10{,}000$ and for HM we take $t=3$
when using $\epsilon=3$ otherwise we use the optimal $t$ values as found in
Figure LABEL:fig:group3_t.

In Figure LABEL:fig:group4_eps, we vary $\epsilon\in\\{0.5,1,\dots,5\\}$,
fixing $d=1024,n=100{,}000$. We can see for low $\epsilon$ values (high
privacy) that the techniques have comparable MSE values. As $\epsilon$
increases we can begin to see a split between the Hadamard technique (HR, HM)
and the local hashing methods, with local hashing methods having lower MSE.
Similarly in Figure LABEL:fig:group4_d, we vary
$d\in\\{2^{2},\dots,2^{11}\\}$. We can see that as $d$ increases the split
between the Hadamard and local hashing methods is more apparent. For very
small domain sizes the techniques are comparable, but as $d$ grows large the
local hashing methods have lower MSE. Finally we look at the total time taken
by the oracles as we vary $d$ in our experiments shown in Figure
LABEL:fig:group4_time. As noted previously, the local hashing methods are slow
as $d$ increases and we can see the comparison in time between OLH and FLH
with $k^{\prime}=10,000$. HR is the fastest technique, followed by HM and then
OUE. We can clearly see that as the domain size increases, the local hashing
methods are much slower.

It should be noted that in these comparisons we use HM with $t=3$ where the
aggregation time takes 3x longer than HM with $t=1$, since each user sends
three Hadamard coefficients. For $t=1$, HM takes almost exactly the same total
time as HR and there is little difference in time between them, except that
HM’s aggregation time increases as $t$ does. Out of all these techniques, OUE
is the only one with $O(d)$ communication cost and as $d$ grows further OUE
would struggle to be competitive with the other techniques.

Figure 6. Varying the regularisation parameter $(\alpha)$ in RAPPOR, fixing
$\epsilon=3,d=100{,}000,n=1{,}000{,}000$

(a)

(b)

(c)

Figure 7. Varying Sketch Vector Size $(c)$ and Number of Hash functions $(r)$
with different sketch methods
Figure 8. MSE (left) and $k$-MSE (right) comparison of Bloom-filter, Sketch
and Count-Sketch with Direct Encoding (DE). Fixing $\epsilon=3,r=16,c=m=128$

###  7.5. Domain Reduction Methods

Next, we experimentally compare domain reduction methods (introduced in
Section 4) which we combine with various FOs to allow us to efficiently
estimate frequencies over very large domains. We devote particular attention
to the following combinations:

  * •

Count-Sketch (CS) and Count-Min Sketch (CM) methods with FLH, Hadamard
Response (HR), HadamardMech (HM)

  * •

Bloom filters with Direct Encoding (DE) (the basis of the RAPPOR system
(Erlingsson et al., 2014))

  * •

Apple’s proposal of Count-Mean Sketch (CMS) as CM (Mean) sketch with OUE and
Hadamard Count Mean Sketch (HCMS) as CM (Mean) with HM ($t=1$) (app, 2017)

In these experiments we now increase the number of users to $n=1{,}000{,}000$
and vary $d$ from $10{,}000$ to $100{,}000$ in increments of $10{,}000$. We
again use synthetic data following a Zipf distribution. Along with reporting
MSE we also consider $k$-MSE which calculates the MSE over a subset of the
top-$k$ most frequent elements in the domain. We fix $k=50$ in all
experiments. We show this metric because as $d$ grows large, there will be
many elements in the domain that have small (or close to $0$) frequency.
Hence, if a frequency oracle was to underestimate frequencies (or predict
close to $0$) it would have a (misleadingly) low MSE but a very high $k$-MSE
as it would fail to accurately estimate the top-$k$ frequent elements. As we
will see, there is often a trade-off in large domains between achieving a low
MSE and a low $k$-MSE.

####  7.5.1. Varying Parameters in RAPPOR

The RAPPOR method combines Bloom filters with the Direct Encoding (DE) FO. One
difficulty in using Bloom filters is that the estimation is not as
straightforward as using sketches. The RAPPOR paper describes an estimation
method using regularised regression for both selecting domain elements that
could be frequent (via LASSO) and then estimating the frequencies of these
elements (via non-negative ridge regression).

In Figure 6, we run an experiment varying the regularisation parameter
$\alpha$ used in the ridge regression that estimates frequencies from the
Bloom filters. We fix $n=1{,}000{,}000,\epsilon=3,d=100{,}000$ and plot three
lines, each of which corresponds to a Bloom filter of size $m=32,64,128$. We
can immediately note two things: Increasing the size of the Bloom filters
results in lower MSE and that the correct choice of regularisation can lead to
improved MSE performance. The best regularisation value found was $0.005$ but
values in the range $0.005-0.01$ have very similar MSE values. As we increase
or decrease the regularisation from this point the MSE of RAPPOR begins to
suffer. The optimal value of regularisation is likely to be specific to the
distribution of the data that is being collected, which may not be known in
advance.

####  7.5.2. Varying Sketch Parameters in Sketch Oracles

One of the most popular deployments of sketch-based LDP frequency oracles are
Apple’s (H)CMS methods (app, 2017). The CMS technique can be thought of as a
frequency oracle that combines a Count-Mean sketch with the OUE frequency
oracle while its Hadamard variant is effectively the Hadamard Mechanism (HM)
oracle combined with Count-Mean sketch. As explained in Section 4, sketch
methods consist of $r$ hash functions that map elements to a row of a $r\times
c$ sized array. We therefore first look at the effect of varying these sketch
parameters $r,c$. We consider both Count-Sketch (CS) and Count-Min sketch (CM)
methods with variants that take either the mean, median or minimum of the
sketch matrix to produce final frequency estimates.

In Figure 7 we vary the choice of sketch method and fix the frequency oracle
to be FLH with $k^{\prime}=500$. We run experiments with $n=1{,}000{,}000$
users, $d=100{,}000$ and $\epsilon=3$. Figures LABEL:fig:group5_sketch_m_mse
and LABEL:fig:group5_sketch_m_k_mse show both the MSE and $k$-MSE of the
sketch methods while varying the size of the sketch vector
$c\in\\{32,64,\dots,2048\\}$ and fixing the number of hash functions to
$r=32$. We can see that in general CM methods perform slightly better than CS
methods, with a lower MSE. As we would expect, we see that increasing the size
of the sketch vector $c$ decreases the MSE across all methods. Figure
LABEL:fig:group5_sketch_m_k_mse shows an interesting trend in $k$-MSE for
sketches that take the minimum estimate: while fixing $r$ and increasing $c$,
the $k$-MSE begins to slowly increase and after the $c=256$ point the median
and mean methods begin to outperform the minimum.

In Figure LABEL:fig:group5_sketch_k_mse we repeat the same experiment but vary
the number of hash functions, $r\in\\{2,4,\dots,512\\}$ and fix $c=1024$. We
only present the $k$-MSE graph for varying $r$, since we found that the MSE
graph for varying $r$ showed very similar behaviour to Figure
LABEL:fig:group5_sketch_m_mse with a general decrease in MSE as the parameter
is increased. Broadly, we see a similar pattern to the previous
experiment—with CM methods slightly outperforming the CS methods. However, the
notable difference is in the $k$-MSE graph, where increasing $r$ results in
worse $k$-MSE performance for sketch methods that take the minimum as their
estimate.

Figure 9. MSE (left) and k-MSE (right) comparing Count-Median sketch methods
(with FLH, HR, HM) and Apple’s CMS/HCMS while varying $d$, fixing
$\epsilon=3,n=1{,}000{,}000,r=32,c=1024$

####  7.5.3. Comparing Sketch Estimation Methods

Across all of Figure 7 we can see that the choice of sketch estimate quite
clearly affects both the MSE and $k$-MSE performance. We see that using
sketches that take the Minimum as the estimate have the lowest MSE followed by
Mean and then the Median. However, the picture becomes more complex when we
look at the $k$-MSE plot in Figure LABEL:fig:group5_sketch_k_mse. Here the
results are somewhat reversed, with Min methods performing the worst followed
by Median and then Mean. As mentioned, as the number of hash functions $r$ is
increasing, the $k$-MSE of the Minimum methods begin to increase.

The likely reason for this is that as the number of hash function increase,
the number of users are effectively split into cohorts where their data
contributes to a single row of the sketch matrix. This means as $r$ increases,
the user base is split between $r$ frequency oracles each of which estimate a
single row in the sketch. This will lead to noisy estimates that underestimate
the true frequencies. So by taking the minimum of these rows we underestimate
the true frequencies. Generally, it seems that for small $r$ and large $c$ the
Mean and Median are comparable and achieve the best balance between MSE and
$k$-MSE.

####  7.5.4. Sketching vs Bloom Filters

In Figure 8 we compare the various sketching methods with the Bloom filter
approach (implemented via RAPPOR). Since RAPPOR can be thought of as a Bloom
filter combined with the Direct Encoding (DE) frequency oracle, we also
compare sketching methods that use DE as their oracle. To ensure a fair
comparison we fix the sketch sizes to $r=16,c=128$ to match that of a Bloom
filter with size $m=128$. For the Bloom filter we also use $k=2$ hash
functions and $8$ cohorts. We fix the regularisation parameter to be
$\alpha=0.005$ as informed by our previous experiments. We fix
$\epsilon=3,n=1{,}000{,}000$ and vary $d$. We observe generally that in this
small sketch regime the Bloom filter approach is certainly competitive. It
achieves MSE that is similar CS (Min) and begins to outperform it as $d$
increases. This also occurs when we look at the $k$-MSE.

Interestingly, when we look at this small sketch regime taking the minimum as
the sketch estimate leads to better performing oracles in both MSE and
$k$-MSE, different to what we saw for much larger sketches. Since we have a
very large domain and a small sketch, there is likely to be many collisions in
the sketch. Taking the minimum seems to be more robust to collisions in the
sketch whereas the mean and median struggle with this on smaller sketches.
There are limitations to the Bloom filter approach. The first is that using DE
restricts us to using relatively small Bloom filters since the variance of DE
scales in the size of the domain. The second, as previously referenced, is the
difficult choice of regularisation and the complicated estimation procedure.
Consequently, we drop the Bloom filter approach from subsequent experiments.

Figure 10. MSE while varying $\epsilon$, fixing
$d=100{,}000,n=1{,}000{,}000,r=32,c=1024$ Figure 11. MSE (left) and $k$-MSE
(right) of various post-processing techniques on FLH, varying $d$ and fixing
$\epsilon=3,n=1,000,000$

####  7.5.5. Comparing Sketching Using Different Frequency Oracles

We now compare the sketching methods while varying the frequency oracles used
to privatise data. Since our previous experiments have shown that Count-Sketch
methods perform slightly worse, we omit them from our comparison. Instead we
compare Apple’s CMS and HCMS (as CM (Mean) OUE and CMS (Mean) HM ($t=1$),
respectively) and CM (Median) with FLH, HR and HM as its oracle. For HM we fix
$t=3$ when $\epsilon=3$ otherwise we use optimal $t$ values found from our
previous experiments. We also plot CM (Mean) HM ($t=3$) to see if there is
direct improvement over Apple’s HCMS. The experimental results are shown in
Figure 9 and 10 where we fix $r=32,c=1024,n=1{,}000{,}000$.

In Figure 10 we vary $\epsilon$ while fixing $d=100{,}000$. We can see that
for low $\epsilon$ values the Mean methods outperform the other sketches but
as $\epsilon$ increases CM with FLH has the best performance. Out of the
frequency oracles used for CM, it seems that FLH performs the best followed by
HR and then HM. As we have seen, when we increase $k^{\prime}$ in the FLH
oracle the MSE is reduced. We have chosen $k^{\prime}=500$ to ensure that
sketching with FLH, HR and HM take roughly the same amount of total time.
Increasing $k^{\prime}$ will make the overall oracle slower but could increase
accuracy further. In general, all of the sketching methods we have compared
have almost identical client and server-side processing times.

In Figure 9 we fix $\epsilon=3$ and vary the domain size. We can see that CM
(Mean) OUE and CM (Median) FLH achieve the lowest MSE and when looking at the
$k$-MSE it seems that in general the CM (Median) methods have lower $k$-MSE
than CM (Mean). We also note that CM (Mean) with HM $(t=3)$ outperforms taking
CM (Mean) HM $(t=1)$ and therefore improvements in Apple’s current deployment
of HCMS could be made by modifying their protocol to sample more Hadamard
coefficients. Additionally, the CM sketching oracles have a logarithmic
communication cost of $O(\log c)$ which for large sketches is much smaller
than CM (Mean) OUE’s $O(c)$ communication. It seems that the CM (Median) FLH
method has similar performance to CM OUE, yet achieves smaller communication
cost without loss of accuracy. By expressing Apple’s (H)CMS protocols in our
framework we have highlighted some modifications and variations that lead to
increased performance.

###  7.6. Post-Processing Experiments

We next run experiments using the post-processing methods described in Section
6. Figure 11 shows the results, fixing
$d=100{,}000,\epsilon=3,n=1{,}000{,}000$ and varying the post-processing
methods used on a FLH oracle with $k^{\prime}=500$.

Just like sketching, the different post-processing methods also have a
MSE/$k$-MSE trade-off. As we would intuitively expect, non-negative
normalisation results in better MSE and $k$-MSE performance than performing no
post-processing at all. Additive Normalisation and Threshold Cut post-
processing strike a balance between MSE and $k$-MSE. For additive
normalisation it seems that the combination of rounding negative values and
normalising individual estimates to match the total number of users helps to
decrease MSE while still retaining accurate top-$k$ estimates ensuring a
relatively low $k$-MSE. For thresholding, we achieve a lower MSE by rounding
values after the threshold down to $0$ but during this process we seem to
round down certain top-$k$ values which inflates its $k$-MSE. For large
domains projecting the results onto a probability simplex pushes the estimates
down, resulting in a low MSE but a very high $k$-MSE since frequent domain
items have their estimates reduced.

Overall it seems that additive normalisation has the best balance between
MSE/$k$-MSE, but non-negative normalisation also performs well and is clearly
a better choice than no post-processing. Note that for sketching oracles, we
can consider two areas to perform post-processing, one that can be performed
on the internal frequency oracles that estimate counts in sketch rows, and a
second post-processing step over the entire sketch oracle itself when
estimating the domain. Our tests showed that taking combinations of additive
normalisation in both of these give the best balance of $k$-MSE/MSE, but we
omit these plots here.

Table 2. Top 20 HH + FO Combinations on the AOL dataset based on F1 scores. The standard deviations of F1 scores are presented in parentheses. The maximum values of precision, recall and F1 are in bold.  | Heavy Hitter | Sketch Method | Frequency Oracle | Parameters | Precision | Recall | F1 Score
---|---|---|---|---|---|---|---
1 | PEM | CM (Median) | FLH | T=20, r=1024, c=2048 | 0.822 | 0.74 |  0.779 (0.052)
2 | PEM | CM (Median) | FLH | T=20, r=128, c=1024 | 0.789 | 0.74 | 0.763 (0.051)
3 | PEM | CM (Mean) | OUE | T=10, r=1024, c=2048 | 0.626 | 0.96 | 0.757 (0.057)
4 | PEM | CM (Mean) | OUE | T=20, r=32, c=1024 | 0.782 | 0.72 | 0.749 (0.025)
5 | PEM | CM (Median) | FLH | T=10, r=1024, c=2048 | 0.758 | 0.74 | 0.748 (0.061)
6 | PEM | CM (Mean) | OUE | T=20, r=128, c=1024 | 0.786 | 0.70 | 0.738 (0.035)
7 | PEM | CM (Mean) | OUE | T=10, r=32, c=1024 | 0.782 | 0.70 | 0.738 (0.025)
8 | PEM | CM (Median) | FLH | T=20, r=32, c=1024 | 0.792 | 0.68 | 0.731 (0.036)
9 | PEM | CM (Mean) | OUE | T=20, r=1024, c=2048 | 0.766 | 0.70 | 0.73 (0.029)
10 | PEM | CM (Median) | FLH | T=10, r=32, c=1024 | 0.744 | 0.70 | 0.721 (0.055)
11 | PEM | CM (Median) | FLH | T=10, r=128, c=1024 | 0.552 | 0.88 | 0.677 (0.031)
12 | PEM | CM (Mean) | OUE | T=10, r=128, c=1024 | 0.541 | 0.88 | 0.669 (0.076)
13 | PEM | CM (Mean) | HM (t=1) | T=20, r=128, c=1024 | 0.616 | 0.60 | 0.607 (0.098)
14 | PEM | CM (Mean) | HM (t=1) | T=20, r=1024, c=2048 | 0.580 | 0.56 | 0.569 (0.093)
15 | PEM | CM (Mean) | HM (t=1) | T=10, r=32, c=1024 | 0.551 | 0.54 | 0.545 (0.1)
16 | PEM | CM (Mean) | HM (t=1) | T=20, r=32, c=1024 | 0.516 | 0.48 | 0.496 (0.052)
17 | SFP | CM (Mean) | OUE | T=10, r=128, c=1024 | 0.967 | 0.32 | 0.465 (0.12)
18 | SFP | CM (Mean) | OUE | T=10, r=1024, c=2048 | 0.960 | 0.30 | 0.45 (0.065)
19 | SFP | CM (Mean) | OUE | T=20, r=1024, c=2048 | 1.000 | 0.28 | 0.436 (0.051)
20 | PEM | CM (Mean) | HM (t=1) | T=10, r=128, c=1024 | 0.322 | 0.62 | 0.424 (0.046)

##  8\. Heavy Hitter Experiments

We conclude our experiments by applying the frequency oracles to the heavy
hitter (HH) problem. We consider three popular implementations of LDP HH
techniques, Apple’s Sequence Fragment Puzzle (SFP) (app, 2017), and two from
the Hierarchical Search framework: Prefix Extending Method (PEM) (Wang et al.,
2019) and TreeHistogram (TH) (Bassily et al., 2020; Bassily and Smith, 2015).

Heavy Hitters Dataset. We base our HH experiments on a subset of the 2006 AOL
search query dataset. The AOL dataset, among other things, contains both
search queries by users and (depending on if they clicked on a website) the
URL of the website they chose. We focus only on the clicked URLs and not the
search queries themselves.444The AOL dataset can be accessed here:
http://www.cim.mcgill.ca/~dudek/206/Logs/AOL-user-ct-collection/ This forms a
dataset of $1{,}935{,}614$ URLs, with $383{,}467$ unique URLs. The most
popular among them are websites like google.com, yahoo.com and myspace.com.
The goal of applying HH techniques in this scenario is to discover the most
popular URLs that were clicked by users. Additional pre-processing was
performed to strip website URLs of leading domain prefixes like www. and
https:// but domain suffixes like .com were left untouched.

###  8.1. Experiment Parameters

We run combinations of heavy hitter search protocols (PEM, SFP, TH) with
various sketch-based frequency oracles with the goal of finding the best
combination of frequency oracle and heavy-hitter search technique for this
dataset. We run the HH + FO combinations to attempt to find the top-10 most
frequent URLs in the dataset. We privatise all $1{,}935{,}614$ URLs with a
privacy budget of $\epsilon=3$, which for SFP and TH is split evenly into
$\epsilon=1.5$ for both the fragment and word estimator. We measure the
performance via the precision, recall and the associated F1 score calculated
by comparing the strings found to the true top-10 frequent URLs in the
dataset.

We truncate all URLs in the dataset to a length of $6$ characters. Across all
the HH protocols, we fix the length of prefixes that are used to build up the
HHs to be a length of $2$ characters. We run two variations of the HH
protocols for $T=10,20$ where $T$ denotes the number of top-T prefixes that
are used to build up the heavy hitters at each stage of the protocol. While
each of the three HH protocols differ in how they do this, all of them
estimate some combination of top-$T$ prefixes from which they build up their
HH predictions. For the frequency oracles we choose Apple’s CMS (as CM (Mean)
+ OUE), HCMS (as CM (Mean) + HM) and CM FLH (Median) to compare with current
HH deployments that use (H)CMS. We use three sketch-size variants $(r,c)$:
$(32,1024),(128,1024),(1024,2048)$. Overall, this results in a total of $54$
combinations, each of which is repeated $5$ times with the results averaged.

###  8.2. Results

In Table 2, we present the top-20 combinations of HH + FO sorted by their F1
score. The most striking result is that PEM dominates SFP and TH in terms of
the F1 score, no matter the oracle that we choose. These results corroborate
similar experiments in (Wang et al., 2019), however they choose to just fix
the oracle to be OLH and compare PEM with slight variations of hierarchical
search techniques.

Focusing on the PEM results, CM FLH has the best F1 performance for the two
larger sketch-size variations with $T=20$, followed by CMS (i.e., CM (Mean)
with OUE) with the largest sketch-size variant. While the F1 scores of using
CMS and CM FLH with PEM on this dataset are very comparable, the most notable
result here is that we can achieve slightly better F1 scores than CMS with
both much smaller sketch-sizes and much lower communication cost by using CM
FLH. To achieve an equivalent communication cost we would have to use HCMS
which as we can observe performs notably worse than CM FLH.

For SFP, the table suggests that combining it with CMS outperforms SFP with CM
FLH. Looking back at Figure 10, we can start to see why this would be the
case. Since SFP splits the privacy budget between estimating fragments and
words, then both oracles in our experiments use $\epsilon=1.5$. Our
experiments in Figure 10 show that for lower privacy budgets, CMS outperforms
CM FLH and this is the likely cause to the better performance of SFP + CMS.
Just like for PEM, we see that CM FLH outperforms HCMS when combined with SFP
as well.

##  9\. Concluding Remarks

We have seen that methods for collecting frequency statistics under the model
of local differential privacy can be practical and accurate. Although there
are many variations, by viewing them through our four layer model, it is
possible to navigate the large space of possibilities to identify the most
effective combinations. Our experimental findings can be summarized as
follows:

Frequency oracles. For small domains, up to hundreds or thousands in size,
methods based on local hashing are most effective. Our proposed “fast” local
hashing (FLH) approach is several times faster than optimal local hashing
(OLH), without much loss of accuracy. For larger domains, Hadamard-based
approaches, particularly Hadamard response (HR) is almost as accurate but can
be orders of magnitude faster.

Domain size reduction. For larger domains, domain size reduction via sketches
works well. Count-min sketch based approaches seem marginally preferable to
Count sketch ones, using estimation with median or mean performing best across
a range of domain sizes up to hundreds of thousands.

Post-processing. Post-processing reduces error, with the method of rounding
negative values to zero being simplest, while adding a small normalization
constant proving a good choice when feasible.

Heavy hitters. Combining all these to solve heavy hitters problems, the
combination of sketching with local hashing seems the best overall, though
other similar parameter settings do almost as well. The hierarchical encoding,
in particular the prefix encoding method (PEM) clearly dominates other
choices.

Last, we note that while we focus on frequency statistics, it is natural to
extend these methods to related questions, such as quantiles and range
queries, using transformations that are by now standard (Cormode et al.,
2019).

###### Acknowledgements.

This work was supported by European Research Council grant ERC-2014-CoG 647557
and the UK government.

## References

  * (1)
  * app (2017) 2017\.  Apple differential privacy technical overview.  https://www.apple.com/privacy/docs/Differential Privacy Overview.pdf.
  * Acharya et al. (2018) Jayadev Acharya, Ziteng Sun, and Huanyu Zhang. 2018.  Communication Efficient, Sample Optimal, Linear Time Locally Private Discrete Distribution Estimation.  _CoRR_ abs/1802.04705 (2018).  arXiv:1802.04705 http://arxiv.org/abs/1802.04705
  * Acharya et al. (2019) Jayadev Acharya, Ziteng Sun, and Huanyu Zhang. 2019.  Hadamard Response: Estimating Distributions Privately, Efficiently, and with Little Communication. In _The 22nd International Conference on Artificial Intelligence and Statistics_ _(Proceedings of Machine Learning Research)_ , Vol. 89. PMLR, 1120–1129.  http://proceedings.mlr.press/v89/acharya19a.html
  * Bassily et al. (2020) Raef Bassily, Kobbi Nissim, Uri Stemmer, and Abhradeep Thakurta. 2020.  Practical Locally Private Heavy Hitters.  _J. Mach. Learn. Res._ 21 (2020), 16:1–16:42.  http://jmlr.org/papers/v21/18-786.html
  * Bassily and Smith (2015) Raef Bassily and Adam D. Smith. 2015.  Local, Private, Efficient Protocols for Succinct Histograms. In _Proceedings of the Forty-Seventh Annual ACM on Symposium on Theory of Computing_. ACM, 127–135.  https://doi.org/10.1145/2746539.2746632
  * Bloom (1970) Burton Bloom. 1970\.  Space/time trade-offs in hash coding with allowable errors.  _Commun. ACM_ 13, 7 (July 1970), 422–426.
  * Charikar et al. (2002) M. Charikar, K. Chen, and M. Farach-Colton. 2002.  Finding Frequent Items in Data Streams. In _Procedings of the International Colloquium on Automata, Languages and Programming (ICALP)_.
  * Cormode et al. (2018) Graham Cormode, Somesh Jha, Tejas Kulkarni, Ninghui Li, Divesh Srivastava, and Tianhao Wang. 2018\.  Privacy at Scale: Local Differential Privacy in Practice.  Tutorial at SIGMOD and KDD.
  * Cormode et al. (2019) Graham Cormode, Tejas Kulkarni, and Divesh Srivastava. 2019\.  Answering Range Queries Under Local Differential Privacy.  _Proc. VLDB Endow._ 12, 10 (2019), 1126–1138.  https://doi.org/10.14778/3339490.3339496
  * Cormode and Muthukrishnan (2005) G. Cormode and S. Muthukrishnan. 2005.  An Improved Data Stream Summary: The Count-Min Sketch and its Applications.  _Journal of Algorithms_ 55, 1 (2005), 58–75.
  * Cormode and Yi (2020) Graham Cormode and Ke Yi. 2020.  _Small Summaries for Big Data_.  CUP.  http://cormode.org/ssbd
  * Erlingsson et al. (2014) Úlfar Erlingsson, Vasyl Pihur, and Aleksandra Korolova. 2014\.  RAPPOR: Randomized Aggregatable Privacy-Preserving Ordinal Response. In _Proceedings of the 2014 ACM SIGSAC Conference on Computer and Communications Security_ , Gail-Joon Ahn, Moti Yung, and Ninghui Li (Eds.). ACM, 1054–1067.  https://doi.org/10.1145/2660267.2660348
  * Fanti et al. (2016) Giulia C. Fanti, Vasyl Pihur, and Úlfar Erlingsson. 2016\.  Building a RAPPOR with the Unknown: Privacy-Preserving Learning of Associations and Data Dictionaries.  _Proc. Priv. Enhancing Technol._ 2016, 3 (2016), 41–61.  https://doi.org/10.1515/popets-2016-0015
  * Wang et al. (2018) Ning Wang, Xiaokui Xiao, Yin Yang, Ta Duy Hoang, Hyejin Shin, Junbum Shin, and Ge Yu. 2018.  PrivTrie: Effective Frequent Term Discovery under Local Differential Privacy. In _34th IEEE International Conference on Data Engineering_. IEEE Computer Society, 821–832.  https://doi.org/10.1109/ICDE.2018.00079
  * Wang et al. (2017) Tianhao Wang, Jeremiah Blocki, Ninghui Li, and Somesh Jha. 2017\.  Locally Differentially Private Protocols for Frequency Estimation. In _26th USENIX Security Symposium, USENIX Security_ , Engin Kirda and Thomas Ristenpart (Eds.). USENIX Association, 729–745.  https://www.usenix.org/conference/usenixsecurity17/technical-sessions/presentation/wang-tianhao
  * Wang et al. (2019) T. Wang, N. Li, and S. Jha. 2019.  Locally Differentially Private Heavy Hitter Identification.  _IEEE Transactions on Dependable and Secure Computing_ (2019).  https://doi.org/10.1109/TDSC.2019.2927695
  * Wang et al. (2020) Teng Wang, Xuefeng Zhang, Jingyu Feng, and Xinyu Yang. 2020\.  A Comprehensive Survey on Local Differential Privacy Toward Data Statistics and Analysis in Crowdsensing.  _CoRR_ abs/2010.05253 (2020).  arXiv:2010.05253 https://arxiv.org/abs/2010.05253
  * Wang and Carreira-Perpiñán (2013) Weiran Wang and Miguel Á. Carreira-Perpiñán. 2013\.  Projection onto the probability simplex: An efficient algorithm with a simple proof, and an application.  _CoRR_ abs/1309.1541 (2013).  arXiv:1309.1541 http://arxiv.org/abs/1309.1541
  * Yang et al. (2020) Mengmeng Yang, Lingjuan Lyu, Jun Zhao, Tianqing Zhu, and Kwok-Yan Lam. 2020.  Local Differential Privacy and Its Applications: A Comprehensive Survey.  _CoRR_ abs/2008.03686 (2020).  arXiv:2008.03686 https://arxiv.org/abs/2008.03686