Scalable Secure Computation of Statistical Functions with
Applications to k-Nearest Neighbors
Hayim Shaul, Dan Feldman, Daniela Rus

arXiv:1801.07301v1 [] 22 Jan 2018

January 24, 2018
Abstract
Given a set S of n d-dimensional points, the k-nearest neighbors (KNN) is the problem of quickly
finding k points in S that are nearest to a query point q. The k-nearest neighbors problem has applications
in machine learning for classifications and regression and and also in searching. The secure version of
KNN where either q or S are encrypted, has applications such as providing services over sensitive (such
as medical or localization) data.
In this work we present the first scalable and efficient algorithm for solving KNN with Fully Homomorphic Encryption (FHE) that is realized by a polynomial whose degree is independent of n, the number
of points. We implemented our algorithm in an open source library based on HELib implementation for
the Brakerski-Gentry-Vaikuntanthans FHE scheme, and ran experiments on MIT’s OpenStack cloud.
Our experiments show that given a query point q, we can find the set of 20 nearest points out of more
than a thousand points in less than an hour.
Our result introduces a statistical coreset, which is a data summarization technique that allows
statistical functions, such as moments, to be efficiently and scalably computed. As a central tool, we
design a new coin toss technique which we use to build the coreset. This coin toss technique and
computation of statistical functions may be of independent interest.

1

Introduction

The use of cloud services is increasing dramatically. Indeed, using a service located in the cloud has its
advantages over downloading and installing the service locally. Consider, for example, (i) a map that can
be queried for a site (such as a restaurant or a hotel) near a certain location; (ii) a diagnostics tool for
skin cancer, where a user uploads an image of her skin; or (iii) a recommender system based on previous
activities of similar profiles. Even if the databases are small enough to be installed on the user’s device, a
cloud service can constantly updates its databases and algorithms to have more accurate and more relevant
results, something that is hard for a user to do. The latter two examples, being examples of machine learning
indeed require updates to the databases as more data is collected from willing users, and also updates to
the algorithms as they are improved. In the first example as well, updates are required to the database as
sites are added or removed. One big disadvantage of a cloud service is the loss of privacy for the user. The
user is expected to reveal its data to be able to use the service. In some cases this data may be something
the user does not wish to reveal. For example, the output of a recommender system can also be used to
perform differential pricing on the expense of the user. In other cases, it is the interest of the service not
to know user’s information as it may impose extra regulation on the service. For example, having access to
user medical data may require it to comply with HIPAA regulations.
Multi-party computation (MPC) is used to let two or more parties compute a certain output from their
inputs without disclosing any information (except maybe for the output) to the other party.
Fully Homomorphic Encryption (FHE) is a secret-key-public-key encryption scheme that also supports + and × operations on the ciphertext. To encrypt a message m, the public key is used m0 =
1

Encrypt(m, pub), and decryption is possible only with the secret key m = Decrypt(m0 , sec). In addition we have, Decrypt(Encrypt(x, pub) + Encrypt(y, pub), sec) = x + y and Decrypt(Encrypt(x, pub) ×
Encrypt(y, pub), sec) = x × y. This suggests that any polynomial P (and in fact any algorithm) can be
computed on encrypted messages (x1 , x2 , . . .) by applying + and × operations on encrypted values to get
P(x1 , x2 , . . .). In the context of multi-party computation (MPC), FHE can be used to construct a singleround protocol where one party (the cloud) applies a polynomial on user’s data, such that the cloud does
not learn anything on the user’s data. The result of applying a polynomial P on encrypted input Encrypt(x)
is an encrypted value Encrypt(P(x)) that can only be decrypted by the user with her secret key. Except for
P(x), the user does not learn anything about P.
The implementation of FHE schemes usually involves encoding the data m using some transformation
over a finite ring while adding some noise. Removing the noise and reversing the transformation is possible
only with the secret key. As operations are performed on ciphertexts the noise grows. The user needs to
make sure the generated secret key will be able to remove the maximal amount of noise accumulated during
the computation of P(x). There are two techniques to reduce noise: modulo switching and bootstrapping.
The latter is not CPU efficient, as it involves the processes of decrypting a ciphertext under FHE. Modulo
switching, on the other hand, is very efficient to use, however, the number of times it can be used is set at
the time the key pair is generated. With FHE schemes known to date, as key-pair are generated to support
more modulo-switching, they grow in size, take more time to generate and render operations performed on
ciphertext to have more overhead.
The main challenge when designing an algorithm for FHE is to make it compatible with keys generated
with the least amount of modulo switching. Since the noise accumulated in a × operation is significantly
larger than the noise accumulated in a + operation, the latter is usually ignored and the challenge becomes
to design an algorithm whose equivalent polynomial has a low degree. Although theoretically any algorithm
can be expressed as a polynomial of its inputs, in practice, the naive polynomials are rarely practical.
Hereinafter, computations are assumed to be compatible to the FHE model (arithmetic circuit model),
where we are restricted to use + and × over some finite ring.
Definition 1.1 (Polynomial’s size and depth) A function f : Zsp → Ztp is a polynomial over a finite ring
Zp if it is a multivariate polynomial that maps a vector x of s ring elements in Zp , to a vector y = f (x) of t
ring elements in Zp , which is the result of the polynomial f on the input x. The complexity of the computation
is measured by the degree of the polynomial f , denoted by deg(f ), and the number of multiplications needed
to evaluate it, that is denoted by size, denoted by size(f ).
√
For a general polynomial, P, with deg(P) = n, we have size(P) = n. See [31] for details. However,
for some polynomials,
sparse, smaller bounds exist. For example, over Zp , where p is prime,
 i necessarilyp−1
P p−1not
the polynomial
x
=
(x
+
1)
, whose value is 0 for x = p − 1 and 1 otherwise, has size O(log p).
i
Motivated by FHE problems, we measure the depth and size of the polynomials in our solution. We focus
on reducing their degree and their size.
Comparing encrypted numbers. Unlike comparison model, implementing a comparison in FHE is not
trivial. Simply put, FHE does not provide a comparison operator. Having the ability to compare two
encrypted numbers and getting an unencrypted result would essentially break the security of the encryption
as it leaks information on the ciphertext. However, for two encrypted numbers JxK and JyK it is possible to
construct a polynomial P (x, y) whose value is encrypted and equals 1 if x < y and 0 otherwise.
k-nearest neighbors (KNN) is a problem where given a finite set S ⊂ Rd , with |S| = n and a point
q ∈ Rd we are looking for a subset S 0 ⊂ S, of the k nearest points to q, i.e. |S 0 | = k and maxx∈S 0 kx − qk ≤
minx∈S\S 0 kx − qk. In machine learning the KNN algorithm is used as a classifier, when each point x ∈ S
is assigned a class and the class of q is taken to be the majority of the classes in S 0 . Specifically, KNN can
address the examples mentioned above. Setting S ⊂ R2 to be a set of “sites” and q ∈ R2 to be users location,
KNN is the subset of k nearest sites to q. Setting S ⊂ Rd to be a training set for skin cancer, labeled as
2

“benign” or “malignant”, and q ∈ Rd be a sample to be classified, KNN is a subset of k nearest training
sample to q from which q can be classified as “benign” or “malignant”. Lastly, setting S ⊂ Rd to be profiles,
each paired with a recommendation, and q ∈ Rd to be a user’s profile, KNN is a subset of k most similar
profiles to q which suggests k recommendations to the user.
Despite its simplicity, k-nearest neighbors is a hard problem to solve in FHE. Informally, a polynomial
whose value is the index of one of the k nearest neighbors should somehow consider all data points in S.
Naively, projecting from a comparison-model algorithm, one would expect this to be a polynomial of degree
O(n). Making it harder, is the fact that in the FHE model an equivalent of the < operator needs to be
implemented from + and × operators which naively results in a polynomial with degree O(n · deg(<)), where
deg(<) is the degree of the polynomial evaluating the < operation.
Statistical functions such as the cumulative distribution function (CDF) can be used to solve the KNN
problem. Specifically, for a random variable d = kq − si k, where i is drawn uniformly from {1, . . . , n}, and
CDFd being its CDF then by definition we have, | {si ∈ S | kq − si k < T } | = k, where T = CDFd−1 (k/n).
A closely related set
P ofe functions are moments. Given a data set D = {d1 , . . . dn }, the e-th moment is
di . The first order moment, with e = 1, is the average of D. The variance of D can
defined as µe (D) = n1
be computed from its first two moments V (D) = µ2 (D) − (µ1 (D))2 . For some distributions the CDF can be
computed from the first few moments.
Computing moments in FHE is a hard problem. In fact, even computing the first moment, i.e. the
average, is hard. The main problem is the lack of a division operation in a ring. Even if a multiplicativeinverse exists, multiplying by the inverse does not act as division. For example, computing the average of
−1
2 and 3, Avg(2, 3) = 2+3
we would have (2 + 3) · 2−1 = 6
2 = 2.5, over the ring Z7 . By multiplying by 2
mod 7, where in fact we would like to have (since we are rounding things to integers) Avg(2, 3) =P
2. The
issue of computing an average has been addressed before. In [30] it is suggested to return the pair ( xi , n)
as the output of an average function. This solves the division problem by deferring it to a later step. This
has two obvious disadvantages: (i) it requires an expensive step to perform the division either by an iteration
with aPTuring machine or by embedding a division polynomial and (ii) in most interesting applications
we
P
have
di  max di > Avg(d1 , . . . , dn ) which requires operations on numbers as large as
di which is
inefficient when n is large.
Finding k smallest values in a general unordered array D = {d1 , . . . , dn }, with di ∈ [L], for 1 ≤ i ≤ n, is
the problem of returning mink D, where mink D ⊂ D is a subset such that mink D = k and max mink D ≤
min(D \ mink D). It is believed this is a hard problem to solve practically in FHE. Informally, we expect
the polynomial degree of such a solution to be at least O(n log L) as it at least needs to perform n − 1
comparisons each of O(log L) bits. Setting D1 = D and iteratively computing mi = min1 Di and setting
Di+1 = Di \ {mi }, for i = 1, . . . n can be used to sort D, a problem that has received a lot of attention in
homomorphic encryption [11, 12, 13, 14] and no efficient solution is yet known (see more in Section 2.3) .

1.1

Our Contribution

In this work we show, counter to common belief, that averages and higher moments can be computed
efficiently in FHE and be used efficiently to solve the k-nearest neighbors problem in near real time. More
specifically our contributions are as follows.
First FHE algorithm to compute statistical functions in a practical way by a polynomial whose
degree is O(deg(<)), the degree of a polynomial comparing encrypted numbers, and independent of n. Our
algorithm is also “embarrassingly
parallel”, which makes it easily scalable. The output of our algorithm is
Pn
1
e
an approximation of m
d
.
For
the applications in this paper we consider the cases where m = n and
i=1 i
e ∈ {1, 2}, but it can be extended to any 0 ≤ e ∈ R and m > 0. The output of our algorithm
is a single
P
ciphertext. This has twoPadvantages over previous techniques, whose output was the pair ( dei , m): (i) we
avoid explicitly keeping
dei which is expensive when n is large and (ii) the output ciphertext can be passed
3

to a second step of an algorithm. This is in contrast to prior solutions (e.g. [30]. We expect this technique
will prove efficient in more algorithms.
A novel technique for coin tossing in FHE. Our solution is non-deterministic and relies on “coin
tosses” with probability that depends on an (encrypted) input. To the best of our knowledge, this is the first
efficient implementation of a coin toss. Since coin tossing is a basic primitive for many random algorithms,
we expect our implementation of coin tossing to have a large impact on future research in FHE.
First practical solution to compute k-nearest neighbors in FHE. Our techniques lead to an algorithm and a single-round protocol where the user sends one message whose size is proportional to the
representation of one point q, and receives one message whose size is proportional to the representation of
k points that are nearest to q. The k-nearest neighbors algorithm has application in searching as well as
in machine learning for classification and regression. A practical solution to this problem opens the way to
many practical systems that employ it.
System and experiments results for search. We implemented our algorithms into a system and ran it
on MIT OpenStack cloud. Our experiments, coinciding with our analysis, show that with a single core the
average of 60,000 7-bit ciphertexts can be computed in about two hours. Also, finding 10 nearest neighbors
in a set of 580 hotels in Boston area rounded to a 2D grid of 100 × 100, took about 20 minutes on a single
server.
Open Source Library for Secure-Search with FHE based on HELib [25] is provided for the community [34], to reproduce our experiments, to extend our results for real-world applications, and for practitioners
at industry or academy that wish to use these results for their future papers or products.

1.2

Novel Techniques

Coin tossing in FHE. The notion of coin tossing is somewhat not trivial in arithmetic circuit model, that
can only realize polynomials. Since a polynomial is restricted to + and × operations it is not clear how or why
it is allowed to simulate a coin toss, let alone, when the probability of the toss depends on an encrypted input
value. Unlike traditional RAM model algorithms, polynomials cannot “call” a sub-polynomial and cannot
“choose” a sub-polynomial to embed. Instead, the process of creating such a polynomial should be such
that, independently of the input, a polynomial is constructed. During the construction of the polynomial
randomness is used to decide on the sub-polynomials that are embedded. We describe our algorithms as
pseudo-code that get a random generator R generating random r parameters for coin tossing. This is of course
just for clarity of presentation. Formally speaking, the random numbers are used during the construction
of the polynomial to determine how it is built.
We show an equivalent of a coin toss.
Lemma 1.2 (Coin Toss) Let n be a integer parameter and D a distribution over N with a cumulative
distribution function CDF : N 7→ [0, 1]. There exists a set of polynomials F = {P1 , . . . , Pn }, such that if
P ← F is drawn randomly and x < n then P(x) = 1 with probability CDF (x), and Pi has degree O(log n)
and size O(log n), for every i = 0, . . . , n.
Using coresets to compute statistical functions.
Given n positive integers, x1 , . . . , xn and m, e > 0
Pn
1
e
x
in
an
efficient way that avoids the need of division and
we show how to compute the mean s∗ = m
i=1 i
does not generate unnecessary large numbers. Some interesting cases for our technique are:
Pn
• take m = n, e = 1 to compute the average n1 i=1 xi .
Pn
• take m = n, e = 2 to compute n1 i=1 x2i .

4

We make use of our new coin tossing technique and compute the polynomial s =
(
xe
1 with probability mi
ai =
0 otherwise.

P

ai , where

∗

It is easy to see that E(s) = s∗ . We give a bound on the probability of error P r(|sP− s∗ | > εs∗ ) < 2exp(− εs3 ).
1
xei , enjoys the best of both
Unlike previous technique, e.g. [30], ourPproposed solution for computing m
1
e
worlds: keeping the ring size small p = m
x(i) + O(1), while still having a polynomial with small degree.

1.3

Definitions

Notation.
For a given integer p ≥ 1 we define [p] = {1, · · · , p}. We denote by JxK the encrypted value
of x, which can be a result of encrypting x by calling encrypt(x), or a result of secure computation, for
example by the cloud, in which case x is not known. We also denote by Jx(bits) K the encryption of x given
in binary representation, where each bit is encoded in a different ciphertext. As common in cryptography,
the arithmetic operations in this paper are applied over a finite set of integers in {0, · · · , p − 1}, where the
modulo operator is used to keep the outcome of each operation in this set. Such a set is formally called the
Zp ring.
Definition 1.3 (Ring) The ring Zp is the set {0, ·· · , p − 1} equipped with multiplication
(·) and addition

(+) operations modulo p, i.e., a · b = (a · b) mod p and a + b = (a + b) mod p for every a, b ∈ Zp .
We use the operator a ← [n] for a random variable a and an integer n, to indicate that a is being drawn
$

randomly from [n]. If a is drawn with equal probabilities we use the notation a ←
− [n].

2

Related Work

It has been shown [23, 24, 38] that two parties (Alice and Bob) can compute any function over their secret
inputs over an interactive protocol that reveals no information beyond the function’s output. The major
disadvantage of an interactive protocol is the communication complexity that grows polynomially in the
time to compute the function. Using an FHE scheme such as [6, 22], a single-round protocol can be devised,
where Alice sends Bob her input encrypted with her public key, Bob computes the function using FHE, and
returns the result to Alice. Since Alice’s secret key is needed to decipher the input and output, Bob cannot
learn anything on the input or the output. On the other side, Alice does not learn anything on Bob’s input
except for the output of the function.

2.1

Searching

Prior work on low degree polynomial solving Private Information Retrieval [7, 10, 17, 20, 21, 33] can be used
to solve the k-nearest neighbors by iteratively applying them to to search for potential neighbors in increasing
order. However, these works assume there is at most one matching record. A recent work [2] removed this
assumption and supports any number of matches. In this case, only the first match is returned. This work
also supports a filter that can be given by Alice as part of the (encrypted) input such as t1 < kq − si k < t2
(i.e., points whose distance from q is at least t1 and at most t2 ). Using this filter can reduce the number of
iterations. A slightly modified filter (setting a lower limit on the index of the matched record) can be used to
iteratively get all matches and not only the first one. Nevertheless, this still requires at least k iterations to
report the k nearest neighbors to q. Another recent work [1] assumes an upper bound s ≥ 0 on the number
of matches to an (encrypted) filter and returns all matches in an output whose size is O(s). This work draws
a link between FHE and sketches and uses a sketch matrix to efficiently compact the vector of matches to a
vector of dimensionality O(s). Although, it is possible to return all k records matching a filter kq − si k < t2
in a single round, Alice is still left with the problem of finding the correct value t2 .

5

2.2

Classification

The problem of classification is to preprocess a set of training classified samples so that given a new sample,
its address can be predicted. The problem is split into two steps. The first is the training of a model (either
supervised or not). The second is to determine the classification of a new sample given a trained model. The
k-nearest neighbors problem has been used for classifications in a non-secure settings. For example, in [3, 19]
kNN was used to classify tissues and cancer. For some models the training step has been implemented in
the secure settings. For example, decision trees [29], Naive Bayes [35] and linear discriminant classifiers [18].
Given a trained model, it was shown in [4] how to use FHE to compute medical prediction functions from
patient’s encrypted data. In another work [5] secure protocols were shown to compute classification for
hyperplane decision, Nave Bayes, and decision trees. However, their work uses interactive protocols (i.e.
more than one round).

2.3

Sorting

If we consider the vector of distances D = (d1 , . . . , dn ), where di is the distance from a query point q to a
point si , then k-nearest neighbors is closely related to sorting. Specifically, after sorting the elements of D,
the k smallest values correspond to the k nearest neighbors. On the other hand, an array (or positive values)
can be sorted if we iteratively find the smallest number by setting k = 1 in a 1-dimensional KN N and
remove it from the array. The problem of sorting an array has been greatly studied. for example in [13, 14].
In [11, 12], Cetin et al. surveyed several sorting algorithms implemented in FHE, and compared them to
their arithmetic circuit that sorts n L-bit numbers, given in their binary representation. The multiplicative
depth of their circuit s O(log n + log L), i.e. the underlying polynomial is of degree O(nL). Although their
algorithm was faster than previous result, it still took almost one minute to sort 64 8-bit numbers.

3

Overview

Generic Polynomials We note that for any function f : [p] → [p] there exists a polynomial Pf,p : Zp →
Pp−1
i
Zp such that Pf,p (x) = f (x), for any x ∈ [p]. Moreover, the polynomial Pf,p (x) =
i=0 ai x can be
constructively built by interpolating over all pairs (x, f (x)). We can use Fermat’s little theorem to keep the
degree low, deg(Pf,p ) ≤ p, since xp = x. When the ring size p is clear from the context we omit it and simply
√
write Pf (x). A general p-degree polynomial has size O( p). See [31]. The generalization of this, a k-variate
polynomial Pf,p (x1 , . . . , xk ) has degree (p − 1)k, and size O(pk/2 ).
√
Throughout the paper we use some special polynomials: Psqrt,p (x) = b xc, and Pbiti ,p (x) the i-th bit in
the binary representation of x, i.e.
(
1 The i-th bit in the binary representation of x is set
Pbiti (x) =
0 otherwise.
Coreset is a data summarization C of a set P of items such that f (C, q) approximates f (P, q) for a desired
function f and a query q from a (usually infinite) set of queries. Coreset is a paradigm in the sense that its
exact definition, structure and properties change from paper to paper. Many coresets were recently suggested
to solve main problems in computational geometry (e.g.[9, 15, 32]), machine learning (e.g. [8, 27]), numerical
algebra [37] and graph theory (e.g. [16]). Recently, coresets were used to solve searching problems in FHE
[1, 2].
We make two extensions to coresets. The first extension is we allow the algorithm applied on C be
different than A. I.e. we require that |A0 (C) − A(S)| < εA(S), for some A0 . The second extension, which is
more related to FHE is we are not interested only in cases where |C| < |S|, but also where A0 (C) is easier
to compute in FHE.

6

Definition 3.1 (FHE coreset) Given a vector S = (x1 , . . . , xn ) and an algorithm A, a (δ, ε, A0 , p) coreset
is a vector C = (c1 , . . . , cm ), where |A0 (C) − A(S)| < εA(S) and cj = Pj (C), where Pj is a polynomial of
degree δ, with probability p.
Sketch Matrix We use a sketch matrix M that enables the client to quickly extract χ from the small
coreset (sketch) vector M · χ that is computed on the server. The coreset is based on what is known as a
sketch or disjunct matrix compressing sparse vectors; see [1, 28] and references therein.
k×n

Definition 3.2 (Sketch matrix [1]) Let s, k, n ≥ 1 be integers. A matrix M ∈ {0, 1}
is called an
(s, n)-sketch matrix, if the following holds. There is an algorithm DECODE that, for every vector χ0 ∈ Rk ,
n
returns a binary vector χ = DECODE(χ0 ) ∈ {0, 1} if and only if χ0 = M · χ. The vector χ0 ∈ Rk is called
the s-coreset of the vector χ.
Such a sketch matrix M can be efficiently constructed (offline, independently of the input) using the
result of Indyk, Ngo and Rudra [28]. They suggest a probabilistic and deterministic construction that can
also be decoded in time that is only poly-logarithmic in the length n of the input vector. For simplicity, we
use the following less efficient but deterministic version.
Theorem 3.3 (Corollary C.3 in [28]) Let s, n ≥ 1 be integers. An (s, n)-sketch report matrix M ∈
k×n
{0, 1}
(see Definition3.2) can be computed deterministically in nO(s log s/ log n) time, and k O(1) space,
where k = O(s2 log n). The corresponding decoding algorithm computes DECODE(y) in k O(1) time.

4

Problem Statement

k Nearest Neighbors Given a set of points S = {s1 , . . . , sn } ⊂ Rd , a query point q ∈ Rd and a parameter
0 < k < n, the search problem is to find a subset S ∗ ⊂ S such that |S ∗ | = k and maxσ∗ ∈S ∗ kσ ∗ − qk ≤
minσ∈S\S ∗ kσ − qk.
In the secure version, we assume two parties Alice and Bob, where Bob holds the set S of size n (for
example, the set of all restaurants in the world), and Alice wants to get from Bob a short list (of size k  n)
of restaurants nearest to her location, without having to download the more than O(k) bits from Bob, or
revealing to Bob where she is. In our model, we assume Bob is semi-honest (i.e. he follows the protocol) and
our goal is to prevent Bob from learning anything on the location of Alice. Furthermore, we are interested in
a non-interactive protocol, meaning that Alice sends a single message of size O(1) to Bob and gets a single
reply from Bob of size O(k).

5

k-nearest neighbors Algorithm

In this section we show how to solve the secure k-nearest neighbors problem in an efficient scalable way. We
build upon an algorithm, which is interesting in its own and described in Section 6, to securely compute
statistical functions. We also give three applications that use k-nearest neighbors.

5.1

k Nearest Neighbors (KNN) with Applications

We first start by showing 3 applications of KNN, then we show how the KNN problem can be solved by
finding the k smallest values in an array. Then we provide some applications.
Find k sites near a location Given n sites S = {s(1), . . . , s(n)} ⊂ [L]2 with their associated information
aux(1), . . . , aux(n) ∈ {0, . . . , 2α }, and an encrypted query JqK ∈ [L]2 a call to KNearestNeighborsp,n,k,ε,R,s,aux (JqK, Jq (bits) K)
returns a sketch of the information of the k-nearest sites in S to q.

7

Algorithm 1: KNearestNeighborsp,n,k,ε,R,d,s,α,aux (JqK, Jq (bits) K)
Parameters:
Parameters:
Parameters:
Parameters:
Parameters:
Input:
Input:
Output:

1
2

3
4

A prime p, and integers n  k > 0, and d, α ≥ 1.
A real number ε > 0.
A family of uniform random generator functions R = {Rr : N → [r]}.
√
p
An n-tuple s, where each tuple is a d-dimensional point, where 0 ≤ sδ (i) ≤ d .
An n-tuple aux of auxiliary data associated with s, where each tuple is an α-bit number.
√
A d-dimensional vector JqK = (Jq1 K, . . . , Jqd K), where 0 ≤ qδ ≤ p/d,
(bits)
(bits)
(bits)
The bit representation Jq (bits) K = (Jq1
K, . . . , Jqd
K), where Jqδ
K are the bits of qδ .
A sketch of size O((1 + ε)k) from which the auxiliary data of at least O((1 − ε)k), nearest neighbors cab

/* Denote s(i) = (s1 (i), . . . , sd (i))
*/
√
/* Also denote sδ = (sδ (1), . . . , sδ (n)) an n-dimensional vector, where 0 ≤ sδ (i) ≤ p/d. */
/* Compute an n-dimensional vector of L1 -distances from q to s(i)
*/
JdistancesK := J(0, . . . , 0)K
for each i ∈ 1, . . . , n do
P
(bits)
(bits)
/* Define: kJaK − bkL1 = δ (Jaδ K − bδ + 2(bδ − Jaδ K)isSmallerp (Jbδ
K, Jaδ
K))
*/
Jdistances(i)K := kJqK − s(i)kL1
Jdistances(bits) (i)K := Convert Jdistances(i)K to bit representations using the polynomials Pbiti ,p
/* KSmallestValues computes an indicator vector for the O(k) least distances.
Algorithm 2.

5
6
7
8

9
10
11
12

See

/* Create a sparse return vector using the sparse indicator vector
JχK := KSmallestValuesp,n,k,R (Jdistances(bits) K)
for each i ∈ 1, . . . , n do
for each δ ∈ 1, . . . , α do
(bits)
(bits)
JKNNδ
(i)K := auxδ
(i)χ(i)
/* Use a sketch to return the sparse result vector. See Section 3 and [1] for
details
JResult(bits) K := J((0, . . . , 0), . . . , (0, . . . , 0))K // a (1 + ε)k-dimensional vector of d-tuples
for each δ ∈ 1, . . . , α do
(bits)
JResultδ
K := M(1+ε)×n · JKNNδ K
Return JResultK

8

*/
*/

*/

Classification System Given n samples S = {s(1), . . . , s(n)} ⊂ [L]d with their associated classes aux(1), . . . , aux(n) ∈
{0, 1}, and an encrypted query JqK ∈ [L]d a call to KNearestNeighborsp,n,k,ε,R,s,aux (JqK, Jq (bits) K) returns a
sketch of the classification of the k-nearest samples in S to q. The classification of q is then taken to be the
majority of classifications.
Recommender system Given n profiles S = {s(1), . . . , s(n)} ⊂ [L]d with recommendations associated to them aux(1), . . . , aux(n) ∈ {0, . . . , 2α }, and an encrypted query profile JqK ∈ [L]d , a call to
KNearestNeighborsp,n,k,ε,R,s,aux (JqK, Jq (bits) K) returns a sketch of the recommendations of the k-nearest profiles in S to q. The set of recommendations to q is then taken to be the set of recommendations to profiles
similar to q.
Algorithm Overview Algorithm 1 has (as parameters) integers n  k > 0, a prime, p and an n-tuple,
s, where each tuple is a d-tuple (s1 (i), . . . , sd (i)) denoting the position of the i-th site in d-space, where
√
0 ≤ sx (i), sy (i) ≤ p/d are integers. The algorithm gets as input an encrypted query q = (q1 , . . . , qd ), where
√
0 ≤ qi ≤ p/d denotes the query location given in binary representation as well as a single ciphertext, and

it outputs a sketch of size O(k) from which the user can infer s(i∗ ) | |q − s(i∗ )|L1 ∈ mink {|q − s(i)|} , that
is, the set of sites that are nearest to q, where kp1 − p2 kL1 is the distance in L1 norm.
In Line 3 Algorithm 1 computes the distance Jdistances(i)K = kJqK − s(i))kL1 , for each 1 ≤ i ≤ n. The L1
P
P
norm is computed using Algorithm 5, ka − bkL1 = δ |aδ − bδ | = δ (a−b+2(b−a)isSmaller(b(bits) , a(bits) )).
In Line 4 the distance distances(i) is converted to bit representation. This is done by computing dlog pe
polynomials Pbit1 ,p (distances(i)), Pbit2 ,p (distances(i)), . . .), where Pbiti ,p (x) = 1 if the i-th bit is set in the
binary representation of x. See Section 3. In Line
5 Algorithm 1 calls Algorithm 2 to get an encrypted
0
indicator vector JχK, where χ(i) = 1 if s(i) ∈ mink s and χ(i) = 0 otherwise, for some k 0 = O(k). In Line 8
the algorithm creates a set of sparse vectors by multiplying the elements of aux by χ. This guarantees that
only the auxiliary data of the k nearest neighbors appear in KN N . To be generic, we have aux(i), the
auxiliary associated with the i-th site, be α bits, where α is a parameter. We keep the result of multiplying
the bits of aux by χ in α n-dimensional binary vectors. By the way Algorithm KSmallestValues constructs
χ, we get that it (and also KN N ) is (1 + ε)k-sparse. In Line 11 we multiply KN N by a precomputed
constant sketch matrix M to get a (1 + ε)k dimensional sketch of KNN.
√

p

Theorem 5.1 Given n points s(1), . . . , s(n) ∈ [ d ]d , an n-tuple of auxiliary data aux, such that aux(i) is
√
p
an α-tuple of bits associated with s(i) and a query point q ∈ [ d ]d . Denote by D = {d(1), . . . , d(n)} the set
of L distances, d(i) = ks(i) − qkL1 . For any 0 ≤ k ≤ n, Algorithm 1 outputs a compact representation of
n 1
o
aux(i) | d(i) ∈ minO(k) D of size O(k) of the auxiliary data associated with the O(k) nearest neighbors of
√
q. The polynomials that compute the sketches have degree O(log p) and size O( p log p + nd log p + αn).
The proof is given in Section 7.5.

5.2

Finding k smallest values in a set

As mentioned in the introduction, the problem of finding the k smallest values in a set D is a hard problem
in FHE. We relax this problem by assuming the set D has a known distribution and specifically we show here
a solution for sets with Gaussian distribution. Another relaxation we make, given D and k, is to compute
minκ D where κ ≈ k (see a more formal relation below).
Algorithm 2 Overview Algorithm 2, has as parameters two integers n  k > 0 and a prime p, where
n is the size of the input and k is the size of the output (i.e. number of smallest values to report). The
algorithm gets as input an encrypted n-dimensional vector d ∈ [p]n . The algorithm assumes the values
of d have Gaussian distribution, (d(1), . . . , d(n)) ∼ N (µ, σ 2 ) (where µ and σ are unknown). For many
applications this is a valid assumption. For example, our experiments (see Section 8.1.3) show that the

9

Algorithm 2: KSmallestValuesp,n,k,R (Jd(bits) K)

1

2

3

4
5

6
7
8
9

Parameters: Two integers, n  k > 0, and a prime p.
Parameters: A random generator functions R = {R : N → [n]}.
√
Input:
Encrypted n-dimensional vector Jdbits K, where d ∈ [ p]n .
Output:
An n-dimensional indicator vector JχK, where χ(i) = 1 if d(i) ∈ minκ {d(i)} and χ(i) = 0 otherwise,
where κ ≈ k with high probability (see details below).
P
di . See Section 6 for details.
*/
/* Compute an estimation of n1
JXK := ProbabilisticAveragep,n,n,1,√p,R (JdK)
P 2
/* Compute an estimation of n1
di
*/
JX 2 K := ProbabilisticAveragep,n,n,2,√p,R (Jd(bits) K)
/* Compute estimation of the deviation
*/
Jσ ∗ K := Psqrt (JXK2 − JX 2 K)
/* Compute a threshold for the report
*/
JT ∗ K := JXK + Φ( nk )Jσ ∗ K
J(T ∗ )(bits) K := Compute bit representation of JT ∗ K using Pbiti
/* Compute an indicator vector for the report
*/
χ := (0, . . . , 0) // an n dimensional vector
for i ∈ 1, . . . , n do
Jχ(i)K := SmallerT hanp (J(T ∗ )(bits) K, Jd(i)K)
return JχK

vector of distances indeed has Gaussian distribution. The algorithm returns an (encrypted) indicator vector

n
χ ∈ {0, 1} , such that χ(i) = 1 iff d(i) ∈ minκ {d(1), . . . , d(n)}, where P r |κ − k| < εcn σµ + Φ−1 nk
<



<τ }|
|{i| d(i)−µ
µ4
σ
2
2exp −εO σ + σ2 +µ2 + µ , where Φ(τ ) =
is the cumulative distribution function. We
n
distinct parameters whose values are known to the algorithm from input that are encrypted. Given the
values of µ and σ it easy to compute T . Unfortunately, as discussed in Section 1 these values are hard to
compute in arithmetic circuit model.
P
d(i) by calling ProbabilisticAverage,
In Line 1 the algorithm computes (an approximation to) the average x = n1
√
with e = 1, m = n, l = p (see Section 6 for more details). In Line 2 the algorithm computes the average of
P
√
d(i)2 by calling ProbabilisticAverage, with e = 2, m = n, l = p (see Section 6 for more
squares x2 = n1
q
√
details). In line 3 we approximate the standard deviation by computing JσK = JxK2 − Jx2 K, where · is
n
computed by applying Psqrt (·). In Lines 7-8 the algorithm computes an indicator vector JχK ∈ {0, 1} , where
−1 k
χ(i) = 1 iff Jd(i)K < JxK + Φ ( n )JσK.
p
Theorem 5.2 Given an array D = (d(1), . . . , d(n)), of n elements, where 0 < d(i) < (p), with Gaussian
distribution, D ∼ N (µ, σ 2 ) and a parameter k  n, then KSmallestValues evaluates an indicator vector
χ = (χ(1), . . . , χ(n)), such that χ(i) = 1 if d(i) < maxκ D, where


 



µ
k
µ4
εn
+ Φ−1
< 2exp −εO σ 2 + 2
P r |κ − k| < √
+
µ
,
n
σ + µ2
2π σ
x −µ
|{i| iσ <τ }|
is the cumulative distribution function of standard Gaussian distribution.Also, χ
where Φ(τ ) =
n
can be computed by a polynomial of degree O(log p) and size O(n log p).

The proof is given in Section 7.4.
In the next section we describe how x and x2 are computed efficiently in arithmetic circuit model.

10

6

Algorithms for Computing Moments and Statistical Functions

P
P 2
xi and n1
xi and other moments can be
In this section we show how statistical functions such as n1
computed efficiently under FHE. We believe computing these functions is of independent interest and it has
a broad set of applications expanding beyond k-nearest neighbors.
Algorithm 3: ProbabilisticAveragep,n,m,e,L,R (Jx(bits) K)
Parameters:
Parameters:
Input:
Output:

Integers, L, p, n, m > 0, such that L < p.
e
A uniform random number generator R : N → [d Lm em].
A vector JxK ∈ [L]n , where values are given in binary form using dlog Le bits.
P
2
Le
An encrypted number Jx∗ K such that P r(|1/m x(i)e − x∗ | > ) < 2e−2d m en

e

1
2
3

4
5
6

c := d Lm e
Set JaK := (J0K, . . . , J0K) to be an cn dimensional encrypted vector
for i ∈ 1, . . . , cn do
be
. See Section 6.1
/* P r(CoinT ossL,R(i)1/e (b) = 1) = cm
(bits)
Ja(i)K = CoinTossL,R(i)1/e (Jx
(di/ce)K)
Pcn
Jx∗ K := i=1 Ja(i)K
return Jx∗ K

*/

xi
Algorithm 3 Overview In Line 1 Algorithm 3 computes a resampling constant to ensure cm
< 1 for any
i. See details below. In Line 4 the algorithm uses cn values r1 , . . . , rn drawn uniformly from [cm], which are
used to set
(
1 if (ri )1/e < xbi/cc
ai =
0 otherwise.
P
P e
1
The algorithm then returns
ai as an estimation to m
xi .
n
Theorem 6.1 Let (x(1),
P. . .e, x(n)) ∈ [L] be an n dimensional vector. Also let 0 < e ∈∗ R and m be an
1
integer. Denote µ = m
xi , then Algorithm 3 describes a polynomial whose value is X , that has degree
e
O(log L) and size O(d Lm en log L) such that

P r (|X ∗ − µ| > εµ) < 2exp(−

µε2
).
3

The proof is given in Section 7.3.
Connection to coresets Our algorithm
functions
P e for computing statistical
P is in fact an FHE-coreset
xi , for S = (x1 , . . . , xn ), A0 (C) =
ci , for C = (c1 , . . . , cn ), and
(see Definition 3.1) with A(S) = n1
cj = CoinT ossL (xj ), which is a polynomial of degree L, and |A0 (C) − A(S)| < εA(S), with high probability.

6.1

Coin Toss

CoinToss Overview The formal way to look at CoinTossn,r is as an infinite set of algorithms, with one
algorithm for every pair of (n, r) for any n, r ∈ Z. Setting n and r, the algorithm is uniquely determined
and it gets an encrypted input Jx(bits) K, where 0 < x < n, given by its binary representation, i.e. as a vector
of bits, where Jx(i)K is the i-th bit of x. Since n and r are “embedded” into the algorithm their values are
known (not encrypted) and therefore can be used in a comparison operation. In Line 2 the algorithm tests
whether r > n ≥ x, and returns J0K if it is. Although a polynomial does not have a comparison operator,
11

Algorithm 4: CoinTossn,r (Jx(bits) K)

1
2
3
4

D

Parameters: Two integers r ≥ 0 and n > 0, where r ←− [n] is drawn from [n] with distribution D.
Input:
An encrypted number JxK, such that 0 ≤ x < n is given in its binary form and x(i) denotes its i-th bit.
Output:
Returns a bit JbK, such that P r(b = 1) = CDFD (x).
if r > n then
return J0K
(bits)
Jy
K := encrypt the binary representation of r
return isSmallern (Jx(bits) K, Jy (bits) K)

this comparison is in fact not encoded in the polynomial. The right way to read this line would be: for a
polynomial identified by parameters r > n the polynomial is degenerated and outputs J0K.
In Line 3 the algorithm extract the bits of r and puts them in binary vector y. Again, this is possible
since r is given as a parameter and should be thought of as: for a polynomial with parameter r, we set
constants y(1) . . . y(log n) whose values equal to the bits of the binary representation of r. The algorithm
then calls isSmallern (Jx(bits) K, Jy (bits) K).
x
, take D to be the uniform distribution over 0, . . . , m.
For example, to get a coin toss with probability m
D

Lemma 6.2 The algorithm CoinTossm,r (Jx(bits) K), where 0 ≤ x < m and r ←− [m] is drawn from [m] with
distribution D, returns a bit b, such that P r(b = 1) = CDFD (x), where CDFD (x) = P r(r < x) is the
cumulative distribution function of D, and the polynomial has degree O(log m) and size O(log m).
The proof for this lemma is given in Section 7.2.

6.2

Comparing Encrypted Values

In this subsection we show, for completeness, how to compare two encrypted values, given by their binary
representation.
Algorithm 5: isSmallern (Jx(bits) K, Jy (bits) K)

Parameters: An integer n > 0.
Input:
Two encrypted numbers JxK, JyK, such that 0 ≤ x, y < n are given in their binary form and x(i) and y(i)
Output:
Return J1K iff y < x, and J0K otherwise.

1
2
3
4
5

/* Denote x(i) and y(i) the i-th bit of x and y, respectively
sameP ref ix := (J1K, J0K, . . . , J0K)
for each i = 1, . . . , dlog ne do
JsameBit(i)K := 1 − (Jx(i)K − Jy(i)K)2
Qi
JsameP ref ix(i + 1)K := j=1 sameBit(j)
Pdlog ne
return i=1 (1 − y(i))x(i)JsameP ref ix(i)K

*/

isSmaller Overview The isSmaller algorithm gets a parameter n, and two encrypted input values JxK an
JyK, given by their binary representation, i.e. as a vector of bits, where Jx(i)K and Jy(i)K are the i-th bit of
x and y, respectively. In line 3 the algorithm computes a binary n-dimensional vector sameBit such that
sameBit(i) = 1 iff x(i) = y(i). In line 4 the algorithm computes a binary n-dimensional vector sameP ref ix
such that sameP ref ix(i) = 1 iff at least i − 1 most significant bits of x and y are equal by multiplying
sameBit(j) for j = 1, . . . , i − 1. The first element sameP ref ix(1) is set to J1K. At line 5 the algorithm
Pdlog ne
returns i=1 (1 − y(i))x(i)JsameP ref ix(i)K.
12

Lemma 6.3 The polynomial realized by Algorithm isSmallern (JxK, JyK) equals 1 if x < y and 0 otherwise,
where x and y are given as dlog ne-bit binary number, and it is of degree O(log n), and size O(log n).
The proof is given in Section 7.1.

7

Analysis

In this section we prove the correctness and efficiency of our algorithms. Unlike the algorithms that were
presented top-down, reducing one problem to another simpler problem, we give the proofs bottom up as
analyzing the efficiency of one algorithm builds upon the efficiency of the simpler algorithm.

7.1

Analysis of isSmaller

Lemma 6.3 The polynomial realized by Algorithm isSmallern (JxK, JyK) equals 1 if x < y and 0 otherwise,
where x and y are given as dlog ne-bit binary number, and it is of degree O(log n), and size O(log n).
Proof. Notice that 1 − (x(i) − y(i))2 = 1 if x(i) = y(i) and 0 otherwise. Therefore, sameP ref ix(i) = 1
if x(j) = y(j), for all j < i, effectively meaning x and y have a prefix of at least i bits. Similarly, since
x(i), y(i) ∈ {0, 1}, we have that (1 − y(i))x(i) = 1 only when 0 = y(i) < x(i) = 1.
Let y < x and let i∗ be the index of the first bit that is not equal in x and y, i.e. 0 = y(i∗ ) < x(i∗ ) = 1.
For all i < i∗ we have (1 − y(i))x(i) = 0, since either x(i) = y(i) = 0 or x(i) = y(i) = 1. For all i > i∗ we
have sameP ref ix(i) = 0. Finally, we have (1 − y(i∗ ))x(i∗ )sameP ref ix(i∗) = 1. Therefore isSmallern (x, y)
evaluates to 1.
On the other way, if isSmallern (x, y) 6= 0, then for some i∗ we have (1 − y(i∗ ))x(i∗ ) 6= 0. Since x(i), y(i) ∈
{0, 1} this means x(i∗ ) = 1 and y(i∗ ) = 0 which implies that sameP ref ix(i) = 0, for i > i∗ . Also, we have
sameP ref ix(i∗ ) 6= 0, which means that x(i) = y(i) for i < i∗ . Combining the two, we get that y < x.
The degree of the sameP ref ixn (x, y) polynomial is dominated by the degree of sameP ref ix, which is
O(log n), Similarly, the size is O(log n).
t
u

7.2

Analysis of Coin Toss

We argue that an equivalent of a coin toss that depends on an encrypted input can be achieved.
D

Lemma 6.2 The algorithm CoinTossm,r (Jx(bits) K), where 0 ≤ x < m and r ←− [m] is drawn from [m] with
distribution D, returns a bit b, such that P r(b = 1) = CDFD (x), where CDFD (x) = P r(r < x) is the
cumulative distribution function of D, and the polynomial has degree O(log m) and size O(log m).
Proof. Let F = {P0 , . . . , Pm−1 }, where Pr (x) = isSmallerm (r, x), where isSmallerm (r, x) = 1 if r < x and 0
D
otherwise. Since r ←− [m] is drawn from [m] withP
distribution D, independent to the choice of x, it follows
that P r(CoinTossm,r (JxK) = 1) = P r(r < x) = i<x P r(r = i) = CDFD (x). Following Lemma 6.3, we
have deg(CoinTossm,r ) = O(log m) and size(CoinTossm,r ) = O(log m).
t
u
The proof of the following lemma follows immediately.
Lemma 1.2 (Coin Toss) Let n be a integer parameter and D a distribution over N with a cumulative
distribution function CDF : N 7→ [0, 1]. There exists a set of polynomials F = {P1 , . . . , Pn }, such that if
P ← F is drawn randomly and x < n then P(x) = 1 with probability CDF (x), and Pi has degree O(log n)
and size O(log n), for every i = 0, . . . , n.

13

7.3

Analysis of ProbabilisticAverage

To prove the correctness of Algorithm 3 we first prove this lemma.
n
Lemma
Pn 7.1 Let x ∈ [L] be an n dimensional vector. Also m > max x(i) be a parameter. Denote µ =
1
1 x(i) and set independent Bernoulli random variables:
m
(
1 with probability x(i)
m
ai =
0 otherwise.

Then,

µε2
ai > (1 + ε)µ < exp(−
),
3
X

µε2
Pr
ai < (1 − ε)µ < exp(−
)
2

Pr

X

Pr

X

and


µε2
)
ai − µ > εµ < 2exp(−
3

Proof. Let x ∈ [L]n be an n dimensional vector. Since ai are independent Bernoulli random variables,
P
P
P
2
1
x(i) = µ and by Chernof we have: P r ( ai > (1 + ε)µ) < exp(− µε3 )
we have that E( ai ) = m
P
P
2
and P r ( ai < (1 − ε)µ) < exp(− µε2 ), from which it immediately follows that P r (| ai − µ| > εµ) <
2

2exp(− µε3 ).

t
u

n
Theorem 6.1 Let (x(1),
P. . .e, x(n)) ∈ [L] be an n dimensional vector. Also let 0 < e ∈∗ R and m be an
1
integer. Denote µ = m
xi , then Algorithm 3 describes a polynomial whose value is X , that has degree
e
O(log L) and size O(d Lm en log L) such that

P r (|X ∗ − µ| > εµ) < 2exp(−

µε2
).
3

Proof. Let x ∈ [L]n , m and e be as defined in the theorem. Assume first that max x(i)e < m. If this
e
0
assumption does not hold, set c = d Lm e, define n0 = cn, m0 = cm, x0 ∈ [Le ]n such that x0i = x(bi/cc) and
continue with x0 , m0 and n0 .
Pn e
Pn
x(i)e
1
∗
Denote µ = m
1 xi . Set X =
1 ai , where ai is a Bernoulli variable with P r(ai = 1) = m . By
2

Lemma 7.1 we have that P r (|X ∗ − µ| > εµ) < 2exp(− µε3 ), as required.
We now bound the size and depth of a polynomial that evaluates one Bernoulli random variable with
e
$
e
e
probability P r(ai = 1) = x(i)
− [m] would require
m . The naive approach CoinTossL ,ri (JxK ), with ri ←
e
computing x in binary representation, a thing that increases the degree
√ and size of the polynomial. To
make Algorithm 3 more efficient, we notice that P r(xe < r) = P r(x < e r) and consider the distribution
√
√
√
|{u≤m | d e ue=r0 }|
De,m with P r(r = r0 ) =
, that draws r (non-uniformly) from e m by setting r = e u, where
m
√
$
u←
− [m] is drawn uniformly from {0, . . . , m − 1}. Thus setting ai = CoinTossL,ri (Jx(i)K), with ri = e ui ,
e

$

where ui ←
− [m] is drawn uniformly. We get from lemma 1.2 that P r(CoinTossL,ri (Jx(i)K) = 1) = x(i)
m . Since
Le
∗
X is constructed by adding d m e outputs of CoinTossL,ri the degree of ProbabilisticAverage is O(log L)
e
and its size is O(d Lm e log L).
t
u

7.4

Analysis of KSmallestValues

p
Theorem 5.2 Given an array D = (d(1), . . . , d(n)), of n elements, where 0 < d(i) < (p), with Gaussian
distribution, D ∼ N (µ, σ 2 ) and a parameter k  n, then KSmallestValues evaluates an indicator vector
14

χ = (χ(1), . . . , χ(n)), such that χ(i) = 1 if d(i) < maxκ D, where



 


µ
εn
k
µ4
−1
2
+µ
,
P r |κ − k| < √
+Φ
< 2exp −εO σ + 2
n
σ + µ2
2π σ
x −µ
|{i| iσ <τ }|
where Φ(τ ) =
is the cumulative distribution function of standard Gaussian distribution.Also, χ
n
can be computed by a polynomial of degree O(log p) and size O(n log p).

Proof.
P
P 2
di and µ2 = n1
di be the first two moments and σ 2 = µ2 − µ2 be the variance of D. Also,
Let µ = n1
k
let T = Q( n ), where Q is the quantile function, i.e. | {di | di < T } | = k. For the Gaussian distribution we
have Q( nk ) = µ + Φ−1 ( nk )σ, where Φ is the CDF function of the standard Gaussian distribution, and Φ−1 is
its inverse. Algorithm
2 computes X and X 2 , using Algorithm 3, as estimations to µ and µ2 , respectively.
q
Denote σ ∗ =

2

X2 − X .
02

2

Set ε0 = 2ε + (4ε2 − 6ε) µµ2 , then P r(X 2 > (1 + ε0 )µ2 ) < exp(− µ23ε ). Since, (1 + ε0 )µ2 − (1 − ε)2 µ2 =
(1 + 2ε)µ2 − (1 + 2ε)µ2 , we have:






µ4
µ2 ε02 + 6µε2
2
∗ 2
2
= exp −εO σ +
P r (σ ) > (1 + 2ε)σ < exp −
+µ
.
3
µ2
002

2

Similarly, set ε00 = ε − (ε2 + 3ε) µµ2 , then P r(X 2 < (1 − ε00 )µ2 ) < exp(− µ22ε ). Since, (1 − ε00 )µ2 − (1 +
ε)2 µ2 = (1 − ε)µ2 − (1 − ε)µ2 , we have:






2µ2 ε002 + 3µε2
µ4
∗ 2
2
2
P r (σ ) < (1 − ε)σ < exp −
= exp −εO σ +
+µ
.
6
µ2
Since

√

1 + 2ε < 1 + ε and

002

√

1 − ε < 1 − ε we get that P r(|σ ∗ − σ| < (1 + ε)σ) < exp(− µ2 ε

02

2

+6µε2
)
3

+

exp(− 2µ2 ε 6+3µε ), with ε0 , ε00 as defined above. We therefore have κ−k
= CDF (T ∗ ) − CDF (T ). For
n
∗
−µ
Gaussian distribution we have CDF (T ∗ ) − CDF (T) = Φ(T σ−µ ) − Φ( T −µ
) < Φ( (1+ε)T
) − Φ( T −µ
σ 
σ
σ ). The
inequality being with probability at least 1 − exp −εO σ 2 +

Φ(a + b) < Φ(a) +

√b
2π

µ4
σ 2 +µ2

+µ

. Since Φ0 (x) <

√1
2π

we have

and therefore

Φ

(1 + ε)T − µ
σ




−Φ

T −µ
σ



ε
<√
2π




µ
−1 k
+Φ ( ) .
σ
n

T −µ
k−κ
∗
n = CDF (T ) − CDF (T ), and for Gaussian distribution weget that
 Φ( σ 4 ) −
(1−ε)T −µ
Φ( T −µ
). The inequality being with probability at least 1−exp −εO σ 2 + σ2µ+µ2 +
σ )−Φ(
σ
< √12π we have

Similarly, we have that
Φ( T

∗

−µ
σ )
0

>

Since Φ (x)


Φ

T −µ
σ




−Φ

(1 − ε)T − µ
σ



ε
<√
2π




µ
−1 k
+Φ ( ) .
σ
n

Putting it all together, we get that if we set T ∗ = µ∗ + σ ∗ Φ−1 ( nk ), then for κ = |{d(i) | d(i) < T ∗ }| =
we have:

∗
Φ( T σ−µ )n



εn
P r |κ − k| < √
2π



µ
+ Φ−1
σ

 



k
µ4
2
< 2exp −εO σ + 2
+µ
.
n
σ + µ2
t
u
15


µ .

7.4.1

Extension to other distributions

Algorithm 2 and Theorem 5.2 assume the distribution of X is Gaussian, however, this can be extended to a
large set of distributions. To see what is needed from a distribution let us review where the assumption of
the distribution of X was used.
CDF with bounded derivative. To bound the number of elements in the output of Algorithm 2, we
need some bound on the derivative of the quantile function, |Q0 (x)| < c, for some c ∈ R. More specifically, we
are interested in such a bound in the segment nk − δ < x < nk + δ, where δ depends on ε and the parameters
of the distribution. Informally, this bound means that the error in computing the first two moments of
X does not change by much the number of elements in {di | di < T ∗ } While such a bound exists for all
distributions, the algorithm is impractical for large values of c. For example, for a Bernoulli distribution
with P r(x = 1) = P r(x = 0) = 1/2, we have CDF (0) = 1/2 and CDF (1) = 1, which means the algorithm
can either return about half of X or all of X. Algorithm 1 assumes k  n. For Gaussian distribution we
have Φ0 (x) < √12π , and for nk  21 we even have Φ0 ( nk + δ)  √12π for sufficiently small values of δ.
Distribution parameters easily computed from moments. The second time the algorithm assumes
the distribution is Gaussian N(µ, σ 2 ) is when computing a threshold T using the quantile function T =
Q(k/n). For general Gaussian distribution we have Q(x) = µ + Φ−1 (x)σ where Φ−1 (x) is the inverse of
the standard CDF for Gaussian distribution. To extend to another distribution one needs to show how to
compute the quantile function for that distribution.

7.5

Analysis of KNearestNeighbors
√

p

Theorem 5.1 Given n points s(1), . . . , s(n) ∈ [ d ]d , an n-tuple of auxiliary data aux, such that aux(i) is
√
p
an α-tuple of bits associated with s(i) and a query point q ∈ [ d ]d . Denote by D = {d(1), . . . , d(n)} the set
of L distances, d(i) = ks(i) − qkL1 . For any 0 ≤ k ≤ n, Algorithm 1 outputs a compact representation of
o
n 1
aux(i) | d(i) ∈ minO(k) D of size O(k) of the auxiliary data associated with the O(k) nearest neighbors of
√
q. The polynomials that compute the sketches have degree O(log p) and size O( p log p + nd log p + αn).
Proof. The correctness of Algorithm 1 follows by combining Theorem 5.2 with the compacting method
described in [1]. More specifically, after computing a set of distances D = {d(1), . . . , d(n)}, Theorem 5.2
states that an indicator vector χ of the smallest κ distances can be computed, where κ ≈ k with high
probability (see Theorem 5.2 for exact details). Using the result of [1], a compact vector χ0 ∈ [p]O(κ) can be
created, such that χ can be reconstructed from χ0 , since it is κ-sparse.
Since Theorem 5.2 states that P r(|κ − k| < εc1 n) < 2exp(−εc2 ), where c1 and c2 are constants that
depend on the distribution parameters of D, we get that χ can be reconstructed from χ0 with that probability.
Since computing an L1 -distance function involves d calls to isSmaller and d multiplications, computing
all the distances can be done by a polynomial of degree O(log p) and size O(nd log p). As argued in [1],
since the sketch matrix, M , is a global constant, the multiplication M χ can be done without ciphertext
multiplications at all.
Putting it all together, we get that Algorithm 1 can be realized by a polynomial of degree O(log p) size
√
O( p log p + nd log p + αn).
t
u
Boosting probability To boost the success probability, KNearestNeighbors can be computed b times
independently, thus boosting the success probability. KNearestNeighbors can fail in two ways and the user
can easily detect both failures. Denote c1 to be the constant as defined above. One type of failure is
when KNearestNeighbors returns κ < (1 − εc1 n)k sites, in which case the sketching algorithm decodes the
κ sites correctly. Another type of failure is when KNearestNeighbors returns κ > (1 + εc1 n)k sites, in
which case the sketching algorithm, designed to decode at most (1 + ε)k sites, fails to decode the sites

16

correctly. Of b independent computations of KNearestNeighbors, the probability that all of them failed is
P r(b failures of KNearestNeighbors) = (P r(KNearestNeighbors failed))b .

8
8.1

The System
k-Nearest Sites On A Map

We implemented the algorithms in this paper into a system that finds the nearest hotels to a query location
in Boston. The system was run on MIT’s OpenStack infrastructure. It maintained a database of 580 hotels
in Boston area (given in clear-text). Prior to sending a query to the system a client needs to generate a
pair (P k, Sk) of public key and secret key. To query the system, the client sends an encrypted location JqK
and the public key used for encrypting it. For multiple queries, the client can use the same (P k, Sk) pair
and the public key can be uploaded only once to the system. The encryption of q is done by the client, and
the encrypted location can be sent from a browser running on a smart phone or from a quadcopter. The
system then runs our secure search algorithm on the cloud using the client’s public key to encrypt constants
as needed by the algorithm, and returns an indicator vector of the k = 13 nearest hotels to q. I.e., we set
the auxiliary data to be a single set bit for all hotels. Changing the system to return the name or location
13
taken as the error parameter. The list
of the hosts is easy by changing the auxiliary data. We set ε = 580
(and the query location) can be decrypted only with the private key available only to the user. The server
can update the database of hotels at any time since it is kept in clear-text, for example, by scraping the
Internet or negotiating with the hotels themselves.
8.1.1

Hardware

Our system is generic but in this section we evaluate it on MIT’s OpenStack cloud infrastructure. We used
one of the standard suggested servers, with 12 2.2 GHz Intel Xeon E5-2630 core and 28 GigaByte of RAM.
These cores are common in standard laptops and servers.
8.1.2

Open Software and Security

The algorithms were implemented in C++. HELib library [34] was used for the FHE scheme implementation,
including its usage of SIMD (Single Instruction Multiple Data) technique. The source of our system is open
under the GNU v3 license and can be found in [34]. Our system and all the experiments below use a security
key of 80 bits. This settings can be easily changed by the client.
8.1.3

Experimental Results

In this section we describe the experiments we ran and their results. We ran experiments on a system
implementing Algorithm 1.
Data We tested our system on a list of 580 hotels in the Boston area, given by their GPS coordinates.
As expected the distribution of hotels around Boston is not even. There are many more hotels around
down-town Boston, for example. We drew a grid on the map around Boston and rounded each hotel to the
nearest grid point. We repeated our tests with grid of various sizes: 100 × 100, 80 × 80 and 60 × 60.
The Experiment We ran Algorithm 1, on 580 hotel locations in Boston area, rounding hotel locations
to be on a grid of L × L and taking p = 2L. We drew a random query point, q on the grid and ran
KNearestNeighborsp,n,k,ε,R,d,s,α,aux (q).
For the experiments that had more than 580 points we duplicated the Boston area hotels as many times
as needed. Recall that since the server operates on encrypted data, the running time is oblivious to the
actual coordinates of the query point q or the hotels’ coordinates si , since the server performs the same
operations regardless to the input.
17

Results Our experimental results are summarized in Figure 1. The experiments measured the total time
to: (i) compute the distances from the query point q to each hotel, di = |q − si |; (ii) compute statistical
functions (µ and σ) on di , and (iii) create an indicator vector χ by comparing each di to a threshold computed
from the statistical functions. Each of these steps is “embarrassingly parallel” [36] which means the running
time on each core decreases linearly with the number of cores.
What dominates the slopes of the graphs? Our analysis in Section 7.5 suggests that the time to
compute the k-nearest neighbors is


L log L + n log L
T = O MUL ·
CORES · SIM D
where
• L × L is the grid size to which we rounded the locations of the hotels
• CORES is the number of computation machines (practically, number of core processors that work in
parallel)
• SIM D (Single Instruction Multiple Data) is the amount of integers that are packed into a single
ciphertext. This SIM D factor is a function of the ring size p and the degree of the evaluated polynomial
O(log L); in HELib, this parameter can be read by calling EncryptedArray::size(), see [26].
• M U L is the times for computing a single multiplication operation.
We had CORES = 10 in all experiments, and the other parameters did not change once the grid size was
set. Our formula, therefore, predicts that the graphs should be linear, as they indeed are. The parameters L
and M U L did not change significantly for different grid sizes and the parameter whose value dominated the
slope of the (linear) graph was SIM D = 1288, 1520 and 980, for grid sizes 60 × 60, 80 × 80 and 100 × 100,
respectively. We plotted in dotted lines the linear graph as was predicted by our formula with the various
SIM D values.

Figure 1: Server’s running time (y-axis) on a single machine on MIT’s OpenStack cloud infrastructure, for
different database sizes (x-axis) of KNearestNeighbors (Algorithm 1) on encrypted query. Each colored curve
represents a different grid size. The dotted lines are the graphs that were expected by our formula.

8.2

Computing Statistical Functions

Since we believe our Algorithm 3 for computing
provide experimental results dedicated to this.

1
m

P

18

xei is a useful tool in other applications as well, we

8.2.1

System Overview

For this simple system, we again assumed the user has a pair (P k, Sk) of public and secret key, such that a
message encrypted with P k can be decrypted only with Sk. The input for this system were a list of values
x1 , . . . , xn encrypted with P k by the user. We also assumed
Pn P k is known to the server. The output of the
system was the encrypted value of fe (x1 , . . . , xn ) = n1 i xei which can be decrypted only with the secret
key Sk. Although this alone is not an interesting functionality for a system, we believe that many systems
include a similar computation as a sub system within them. We ran this system as well on MIT’s OpenStack
infrastructure.
P
P 2
xi and n1
xi , i.e. e = 1 (the average) and e = 2 (the
The Experiment We computed the sums n1
average of the squares), the first two moments. We tested with various values of 1, 000 ≤ n ≤ 100, 000, and
ran experiments with xi ∈ [60], xi ∈ [80] and xi ∈ [100]. The server is oblivious to the values of xi since
they are encrypted, and the time of computation is not influenced by the values of xi .
Results Our experimental results are summarized in Figure 2. In Fig. 2Pwe show that indeed the formula
n e
1
correctly predicts the experimental results. Assume we are computing m
1 xi for some values of n, m and
e, where xi < L for some value of L and 1 ≤ i ≤ n, then the overall running time is:


n · log L · M U L
T =O
,
CORES · SIM D
where CORES, SIM D and M U L are as explained in section 8.1.3.

Figure 2: Server’s running time (y-axis)
P on a 1single
P 2 core on MIT’s OpenStack cloud infrastructure, for
computing the first two moments ( n1
xi and n
xi ) (Algorithm 3) of n (x-axis) encrypted values. Solid
lines are for the time to compute the first moment; dashed lines are for the time to compute the second
moment. The dotted lines are the graphs that were expected by our formula. For each moment, the values
were taken from three sets: {0, . . . , 60}, {0, . . . , 80} and {0, . . . , 100}, each represented with a different color.

What dominates the slope of the graph? Although all graphs are linear as expected from our theoretic
analysis, it is interesting to note how their slopes differ. From the formula we gave above, the slope of each
log L·M U L
graph should be proportional to CORES·SIM
D . In our experiments we had CORES = 1 and negligible
differences in log L. Also the differences we measured in the values of M U L were negligible. The values of
SIM D were chosen automatically by HELib upon key generation and varied greatly between experiments
from 192 to 660. To compare our measured results to theory we plotted with dotted lines the linear functions
log L·M U L
1
predicted by our formula: f (x) = xc SIM
D , where c = CORES and SIM D is the SIMD factor that was
set in each experiment. It can be seen that the experiments match our expected formula.
19

References
[1] A. Akavia, D. Feldman, and H. Shaul. Secure database queries on the cloud: Homomorphic encryption
meets coresets, 2017. Submitted.
[2] A. Akavia, D. Feldman, and H. Shaul. Secure search on the cloud via coresets and sketches. arXiv
preprint arXiv:1708.05811, 2017.
[3] A. Ben-Dor, L. Bruhn, N. Friedman, I. Nachman, M. Schummer, and Z. Yakhini. Tissue classification
with gene expression profiles. Journal of Computational Biology, 7(3-4):559–583, July 2004.
[4] J. Bos, K. Lauter, and M. Naehrig. Private predictive analysis on encrypted medical data. 50, 05 2014.
[5] R. Bost, R. A. Popa, S. Tu, and S. Goldwasser. Machine learning classification over encrypted data. In
22nd Annual Network and Distributed System Security Symposium, NDSS 2015, San Diego, California,
USA, February 8-11, 2015, 2015.
[6] Z. Brakerski, C. Gentry, and V. Vaikuntanathan. (leveled) fully homomorphic encryption without
bootstrapping. In Proceedings of the 3rd Innovations in Theoretical Computer Science Conference,
ITCS ’12, pages 309–325, New York, NY, USA, 2012. ACM.
[7] Z. Brakerski and V. Vaikuntanathan. Efficient fully homomorphic encryption from (standard) lwe.
SIAM Journal on Computing, 43(2):831–871, 2014.
[8] V. Braverman, G. Frahling, H. Lang, C. Sohler, and L. F. Yang. Clustering high dimensional dynamic
data streams. arXiv preprint arXiv:1706.03887, 2017.
[9] V. Braverman, A. Meyerson, R. Ostrovsky, A. Roytman, M. Shindler, and B. Tagiku. Streaming kmeans on well-clusterable data. In Proceedings of the twenty-second annual ACM-SIAM symposium on
Discrete Algorithms, pages 26–40. Society for Industrial and Applied Mathematics, 2011.
[10] G. S. Çetin, W. Dai, Y. Doröz, W. J. Martin, and B. Sunar. Blind web search: How far are we from a
privacy preserving search engine? IACR Cryptology ePrint Archive, 2016:801, 2016.
[11] G. S. Çetin, Y. Doröz, B. Sunar, and E. Savas. Depth optimized efficient homomorphic sorting. In
Progress in Cryptology - LATINCRYPT 2015 - 4th International Conference on Cryptology and Information Security in Latin America, Guadalajara, Mexico, August 23-26, 2015, Proceedings, pages 61–80,
2015.
[12] G. S. Çetin, Y. Doröz, B. Sunar, and E. Savas. Low depth circuits for efficient homomorphic sorting.
IACR Cryptology ePrint Archive, 2015:274, 2015.
[13] A. Chatterjee, M. Kaushal, and I. Sengupta. Accelerating Sorting of Fully Homomorphic Encrypted
Data, pages 262–273. Springer International Publishing, Cham, 2013.
[14] A. Chatterjee and I. Sengupta. Searching and sorting of fully homomorphic encrypted data on cloud.
IACR Cryptology ePrint Archive, 2015:981, 2015.
[15] K. L. Clarkson. Coresets, sparse greedy approximation, and the frank-wolfe algorithm. ACM Transactions on Algorithms (TALG), 6(4):63, 2010.
[16] A. Czumaj, C. Lammersen, M. Monemizadeh, and C. Sohler. (1+ ε)-approximation for facility location in data streams. In Proceedings of the twenty-fourth annual ACM-SIAM symposium on Discrete
algorithms, pages 1710–1728. Society for Industrial and Applied Mathematics, 2013.
[17] Y. Doröz, B. Sunar, and G. Hammouri. Bandwidth efficient pir from ntru. In Financial Cryptography
Workshops, pages 195–207, 2014.

20

[18] W. Du, S. Chen, and Y. S. Han. Privacy-preserving multivariate statistical analysis: Linear regression
and classification. In In Proceedings of the 4th SIAM International Conference on Data Mining, pages
222–233, 2004.
[19] S. Dudoit, J. Fridlyand, and T. P. Speed. Comparison of discrimination methods for the classification of
tumors using gene expression data. JOURNAL OF THE AMERICAN STATISTICAL ASSOCIATION,
97(457):77–87, 2002.
[20] Y. Gahi, M. Guennoun, and K. El-Khatib. A secure database system using homomorphic encryption
schemes. arXiv preprint arXiv:1512.03498, 2015.
[21] C. Gentry. A Fully Homomorphic Encryption Scheme. PhD thesis, Stanford University, Stanford, CA,
USA, 2009. AAI3382729.
[22] C. Gentry. Fully homomorphic encryption using ideal lattices. In Proceedings of the Forty-first Annual
ACM Symposium on Theory of Computing, STOC ’09, pages 169–178, New York, NY, USA, 2009.
ACM.
[23] O. Goldreich. Foundations of Cryptography: Volume 2, Basic Applications. Cambridge University Press,
New York, NY, USA, 2004.
[24] O. Goldreich, S. Micali, and A. Wigderson. How to play any mental game. In Proceedings of the
Nineteenth Annual ACM Symposium on Theory of Computing, STOC ’87, pages 218–229, New York,
NY, USA, 1987. ACM.
[25] S. Halevi. Helib - an implementation of homomorphic encryption. https://github.com/shaih/HElib/,
2013.
[26] S. Halevi and V. Shoup. Algorithms in helib. In 34rd Annual International Cryptology Conference,
CRYPTO 2014. Springer Verlag, 2014.
[27] J. Huggins, T. Campbell, and T. Broderick. Coresets for scalable bayesian logistic regression. In
Advances in Neural Information Processing Systems, pages 4080–4088, 2016.
[28] P. Indyk, H. Q. Ngo, and A. Rudra. Efficiently decodable non-adaptive group testing. In Proceedings of
the Twenty-first Annual ACM-SIAM Symposium on Discrete Algorithms, SODA ’10, pages 1126–1142,
Philadelphia, PA, USA, 2010. Society for Industrial and Applied Mathematics.
[29] Y. Lindell and B. Pinkas. Privacy preserving data mining. In Proceedings of the 20th Annual International Cryptology Conference on Advances in Cryptology, CRYPTO ’00, pages 36–54, London, UK,
UK, 2000. Springer-Verlag.
[30] M. Naehrig, K. Lauter, and V. Vaikuntanathan. Can homomorphic encryption be practical? In
Proceedings of the 3rd ACM Workshop on Cloud Computing Security Workshop, CCSW ’11, pages
113–124, New York, NY, USA, 2011. ACM.
[31] M. S. Paterson and L. J. Stockmeyer. On the number of nonscalar multiplications necessary to evaluate
polynomials. SIAM Journal on Computing, 1973.
[32] J. M. Phillips. Coresets and sketches. arXiv preprint arXiv:1601.00617, 2016.
[33] S. S. Roy, F. Vercauteren, J. Vliegen, and I. Verbauwhede. Hardware assisted fully homomorphic
function evaluation and encrypted search. IEEE Transactions on Computers, 2017.
[34] H. Shaul. LiPHE: Library for practical homomorphic encryption., 2017.
[35] J. Vaidya, M. Kantarcöğlu, and C. Clifton. Privacy-preserving naı̈ve bayes classification. The VLDB
Journal, 17(4):879–898, July 2008.
21

[36] B. Wilkinson and M. Allen. Parallel programming: techniques and applications using networked workstations and parallel computers. Prentice-Hall, 1999.
[37] D. P. Woodruff et al. Sketching as a tool for numerical linear algebra. Foundations and Trends R in
Theoretical Computer Science, 10(1–2):1–157, 2014.
[38] A. C.-C. Yao. How to generate and exchange secrets. In Proceedings of the 27th Annual Symposium
on Foundations of Computer Science, SFCS ’86, pages 162–167, Washington, DC, USA, 1986. IEEE
Computer Society.

22

