# Geometric Covering using Random Fields

Felipe Gonçalves  Department of Mathematics, IMPA<EMAIL_ADDRESS>,  Daniel
Keren  Department of Computer Science, University of Haifa
<EMAIL_ADDRESS>,  Amit Shahar  Department of Computer Science,
University of Haifa<EMAIL_ADDRESS>and  Gal Yehuda  Department
of Computer Science, Technion-IIT<EMAIL_ADDRESS>
###### Abstract.

A set of vectors $S\subseteq\mathbb{R}^{d}$ is
$(k_{1},\varepsilon)$-clusterable if there are $k_{1}$ balls of radius
$\varepsilon$ that cover $S$. A set of vectors $S\subseteq\mathbb{R}^{d}$ is
$(k_{2},\delta)$-far from being clusterable if there are at least $k_{2}$
vectors in $S$, with all pairwise distances at least $\delta$. We propose a
probabilistic algorithm to distinguish between these two cases. Our algorithm
reaches a decision by only looking at the extreme values of a _scalar valued_
hash function, defined by a _random field_ , on $S$; hence, it is especially
suitable in distributed and online settings. An important feature of our
method is that the algorithm is oblivious to the number of vectors: in the
online setting, for example, the algorithm stores only a constant number of
scalars, which is independent of the stream length.

We introduce random field hash functions, which are a key ingredient in our
paradigm. Random field hash functions generalize _locality-sensitive hashing_
(LSH). In addition to the LSH requirement that “nearby vectors are hashed to
similar values”, our hash function also guarantees that the “hash values are
(nearly) independent random variables for distant vectors”. We formulate
necessary conditions for the kernels which define the random fields applied to
our problem, as well as a measure of kernel optimality, for which we provide a
bound. Then, we propose a method to construct kernels which approximate the
optimal one.

##  1\. Introduction

The question whether a set of vectors $S$ in Euclidean space is “clusterable”
is of considerable interest in numerous applications. It is also known to be
hard, in its various formulations. For example, the _geometric set cover
problem_ is NP-complete even in dimension two.

We consider the following promise-problem version: given a set
$S\subseteq\mathbb{R}^{d}$, can $S$ be covered by $k_{1}$ balls of radios
$\varepsilon$, or are there at least $k_{2}$ vectors in $S$ with mutual
distances at least $\delta$? The two cases are shown pictorially in Figures
2,2.

We motivate our approach by presenting two basic examples.

###  1.1. Generalized Distinct Count

Consider the following version of the classical Flajolet-Martin (FM hereafter)
algorithm [11], which allows to efficiently, and in constant space, count
distinct items in a data stream. A uniform hash function from the object space
$W$ to the unit interval $h:W\rightarrow I\triangleq[0,1]$ is defined, and the
number of distinct items is estimated by $\frac{1}{\min(h(x))_{x\in S}}$.
Standard “stabilization” techniques (e.g. median of averages) can be applied
to render this estimate more accurate and robust. As the estimate uses only
the minimal value of the hash function, it requires only a constant-size
communication overhead to implement in a distributed setting, as the nodes can
agree on $h$ in advance, and can only broadcast a fixed-size table required to
obtain a stable estimate of the minimum [3].

Assume now that $W$ is not composed of discrete objects, but scalars (real
numbers or vectors with scalar components). In this case, typically, objects
whose corresponding scalars are very close to each other, can be considered
identical, or nearly identical, and it makes little sense to count them as
distinct elements. Instead, we would like to view close vectors as a single
element, and count the tight clusters they form. In this context, clustering
can be seen as a generalization of the Distinct Count problem (“how many
distinct elements are in a given set?”), when generalized to a continuous
space, with “distinct” replaced by “of distance greater than some threshold”.

###  1.2. Inferring on the Size of a Distributed Set

Assume a large set of data vectors are distributed among nodes, and these
nodes wish to determine, with high probability, whether the union of their
data sets can be covered by a ball of radius $\varepsilon$, or whether an
“opposite” condition holds – for example, that the union contains at least two
vectors whose distance from each other is at least $(2+c)\varepsilon$, for
some $c>0$. Further, the decision is to be taken with minimal communication
overhead.

Our proposed algorithm will achieve the desired goal with communication
overhead that is independent of the overall size of the data, and is constant
for each node, regardless of the number of data vectors. Actually, we will
answer a more general question – whether the vectors can be covered by a
certain number of balls, or the “opposite” condition, being the presence of a
set of vectors with pairwise distance greater than a given threshold.

Figure 1.  An example of clusterable data. That is, data that can be covered
with 2 balls with radius $\varepsilon$.

Figure 2. An example of unclusterable data: three vectors with pairwise
distance $\delta>2\varepsilon$.

###  1.3. Our Contribution

We propose a probabilistic approach for the clustering promise problem. Our
main contribution is bringing together ideas from the theory of random fields
and streaming algorithms. More precisely, we show how to use random fields in
order to construct a hash function $h:\mathbb{R}^{d}\rightarrow\mathbb{R}$
with the following properties.

  1. (1)

For every $x\in\mathbb{R}^{d}$, $h(x)$ is a random variable, with identical
distribution.

  2. (2)

Roughly speaking, if $\|x-y\|\leq\varepsilon$, then $|h(x)-h(y)|$ is “small”.
This will ensure that “tight” clusters will be mapped to nearby values. This
property also plays a major role in the _locality sensitive hashing_ (LSH)
scheme [18].

  3. (3)

If $\|x-y\|\geq\delta$ , $h(x)$ and $h(y)$ should be “nearly independent”.
This property ensures that for $x,y$ which are far apart, $h$ behaves like the
hashing used in FM, i.e. it will treat $x,y$ as distinct objects. This
property is not required in LSH, and indeed does not hold for the hash
functions used in its implementation (Section 4).

Intuitively, a hash function satisfying the above conditions can answer the
clustering promise problem: for distant vectors $x,y$, the random variables
$h(x)$ and $h(y)$ are “almost” independent. In this case, similarly to FM
hashing, some simple property of the distribution of $h(x)$ (e.g. its maximum
over some set) carries enough information about the structure of the set (e.g.
its cardinality). On the other hand, for close vectors, the random variables
$h(x)$ and $h(y)$ are highly correlated. This property allows to treat close
vectors as if they are identical.

We show that the function $h$ can be defined using a random field. While we
defer the formal definition for later, a random field over $\mathbb{R}^{d}$
can be thought of as a function that maps each $x\in\mathbb{R}^{d}$ to random
variable $f(x)$ over $\mathbb{R}$. We show explicitly how to construct such
random fields.

There is an inherent tension between Condition 2 and 3: while Condition 2
requires $h(x)$ and $h(y)$ to have “large” covariance if
$\|x-y\|\leq\varepsilon$, Condition 3 requires $h(x)$ and $h(y)$ to have
“small” covariance if $\|x-y\|>\delta$. Ideally, the covariance function of
$h(x)$ and $h(y)$, as a function of the distance, should “drop as fast as
possible” between the distance values $\varepsilon$ and $\delta$. However,
this cannot be achieved by the trivial solution of taking a step function
which equals 1 at $\varepsilon$ and 0 at $\delta$, as it will not define a
distribution, since the corresponding covariance matrix will not be positive
semidefinite; more on this later (Section 6).

We analyze the rate at which the covariance function can decrease, prove lower
bounds, and show how to attain a solution which is near optimal.

###  1.4. Previous Work

A work by Alon et al. [2] investigates the clustering problem in the context
of property testing. There are similarities in the problem definition, in
particular, the goal is to answer whether a set of vector is clusterable, or
“far” from being clusterable. However, there is a key difference between the
property testing setting and the promise problem we consider. The objective in
the property testing setting is to answer the if an “unknown” set
$S\subseteq\mathbb{R}^{d}$ is clusterable or far from being clusterable, while
querying as few vectors from $S$ as possible. More generally, the goal of a
property testing algorithm is testing whether a set $S$ satisfies a property,
or is $\varepsilon$-far from satisfying a property, and the objective is to
make less than $|S|$ queries to $S$ – ideally, the number of queries is
independent of $|S|$ and depends only on $\varepsilon$. The running time, as a
function of $\varepsilon$, is usually less significant. Our goal in this work,
however, is to devise algorithms for answering the clustering promise problem
quickly, while having full access to the dataset.

Additional work in the context of property testing is given in [7]. The
presented approach is also applicable to other symmetric shapes (not just
balls), but it only addresses the problem of covering with a _single_ shape.
We have not been able to find work along the major direction of this paper,
i.e. the application of random fields to the clusterability problem.

For the simpler problem of _distinct count_ (which is equivalent to covering
with balls of zero radii), a well-known approach was presented in FM [11]. One
version of FM is to randomly hash the dataset $S$ into the interval $[0,1]$. A
standard exercise proves that if there are $m$ distinct elements, the
expectation of the minimum value of the hash is $\frac{1}{m+1}$. While our
approach has a similar flavor – we look at the maximal value of a hash
function – it is different in a crucial aspect: FM hashing is oblivious to the
distance between vectors, i.e. even when provided with a very tight cluster,
its output will still depend only on the number of vectors in the cluster;
this is due to the fact that the values of the hash function at any distinct
vector are independent random variables.

The property of mapping nearby vectors to similar values might remind one of
_locality sensitive hashing_ (LSH) method [18], especially the version
presented in [9], which deals with real-valued input. However, while the
requirement “nearby vectors are hashed to similar values” is indeed satisfied
in LSH, the requirement “the hash values are (nearly) independent for far
apart vectors” is not; hence, as will be elaborated later, LSH is
inappropriate for our problem. Still, there appears to be an interesting
relation between LSH and one of the hash functions we propose to use (Section
4).

###  1.5. Paper Structure

After introducing some basic notions in Section 2, we present, in Section 3,
our main contribution and algorithm. Section 4 discusses a new class of random
fields which can also be applied to our problem and that are easier to
construct. Sections 6, 7, 8 are technical in nature, and concern the
construction of optimal kernels, in one and many dimensions. Conclusions are
offered in Section 9.

##  2\. Preliminaries

In this section we fix notations, formally define our problem, and give
standard definitions from the theory of random fields.

###  2.1. Notations

All vector norms are $\ell^{2}$-norms, unless otherwise stated. The inner
product of vectors $x,y$ is denoted $x\cdot y$. For a vector
$v\in\mathbb{R}^{d}$, and $\varepsilon\in\mathbb{R}_{+}$, we denote by
$B(v,\varepsilon)$ a ball of radius $\varepsilon$ and center $v$. In places
where the center is irrelevant, we omit the center from the notation.

###  2.2. Clustering promise problem

We define the clustering promise problem
$\Pi=\Pi_{\text{YES}}\cup\Pi_{\text{NO}}$. Let $S\subseteq\mathbb{R}^{d}$,
$k_{1}\in\mathbb{N}$ and $\varepsilon\in\mathbb{R}_{+}$.

######  Definition 1.

Given a set $S\subseteq\mathbb{R}^{d}$ and parameters $k_{1}\in\mathbb{N}$,
$\varepsilon\in\mathbb{R}_{+}$, we say that $S$ is _$(k_{1},\varepsilon)$
-clusterable_ if there exist $k_{1}$ balls
$\mathcal{B}=\\{B_{1},\ldots,B_{k_{1}}\\}$ of radius $\varepsilon$, such that
for all $x\in S$ there exists a ball $B\in\mathcal{B}$ such that $x\in B$.

We denote by $\Pi_{\text{YES}}=\\{S\subseteq\mathbb{R}^{d}:S\text{ is
}(k_{1},\varepsilon)\text{ clusterable}\\}$. For clarity, we omit the
parameters $\varepsilon,k_{1}$ from the set notation.

######  Definition 2.

Given a set $S\subseteq\mathbb{R}^{d}$ and parameters $k_{2}\in\mathbb{N}$,
$\delta\in\mathbb{R}_{+}$, we say that $S$ is _$(k_{2},\delta)$ -far from
being clusterable_ if there are at least $k_{2}$ vectors in $S$, with all
pairwise distances at least $\delta$.

We denote by $\Pi_{\text{NO}}=\\{S\subseteq\mathbb{R}^{d}:S\text{ is
}(k_{2},\delta)\text{-far from being clusterable}\\}$.

###  2.3. The Relation Between the Conditions $\Pi_{\text{YES}}$ and
$\Pi_{\text{NO}}$

The complementary nature of the conditions is realized by the following Lemma.

######  Lemma 1.

If $S\subset\mathbb{R}^{d}$ has at most $\ell$ vectors with pairwise distance
at least $\delta$, then its minimal cover with balls of radius $\varepsilon$
is of size at most $\ell c(d)\left(\frac{\delta}{\varepsilon}\right)^{d}$,
where $c(d)$ depends only on $d$.

###### Proof.

Suppose there are $\ell$ vectors with mutual distances at least $\delta$ in
$S$, but not $\ell+1$ such vectors. Then, $S$ can be covered by balls of
radius $\delta$ around the $\ell$ vectors (since if these balls do not cover
$S$, there will be a vector at distance at least $\delta$ from all the $\ell$
vectors, a contradiction).

Next, we cover these $\delta$-balls with $\varepsilon$-balls (Fig. 3). The
number of $\varepsilon$-balls required is equal to
$c(d)\left(\frac{\delta}{\varepsilon}\right)^{d}$, where $d$ is the dimension
of the ambient Euclidean space and $c(d)$ a function which depends only on $d$
[26]. Multiplying by $\ell$ completes the proof. ∎

Figure 3. A sketch of the proof of Lemma 1. First, $S$ is covered with balls
of radius $\delta$; then, these balls are covered with balls of radius
$\varepsilon$.

###  2.4. Random Fields

We recall a few definitions relating to the notion of a random field.

######  Definition 3.

(Random Field, [1, 16, 5]). Given a probability space
$(\Omega,\mathcal{F},P)$, a random field $T(x)$ defined on $\mathbb{R}^{d}$ is
a function such that for every $x\in\mathbb{R}^{d}$, $T(x)$ is a random
variable on $(\Omega,\mathcal{F},P)$.

######  Definition 4.

(Covariance of a random field, [1]). The covariance function
$\mathrm{Cov}\left(x,y\right)$ of a random field is defined as

$\mathrm{Cov}(x,y)=\mathop{\mathbb{missing}}{E}[(T(x)-\mathop{\mathbb{missing}}{E}[T(x)])(T(y)-\mathop{\mathbb{missing}}{E}[T(y)])]$

A random field is isotropic and stationary if there exists a function $f$ such
that $\mathrm{Cov}(x,y)=f(\|x-y\|)$ 111In fact, it is possible to define these
two properties separately, but in this work we are only interested in fields
which are both isotropic and stationary.. Since we are interested in
properties which depend only on the distances between the vectors of $S$, we
will only consider such fields.

##  3\. Hash Functions for Geometric Clustering

At the foundation of our algorithm is a random field which, loosely speaking,
maps nearby vectors to similar values, while the images of far away vectors
are mapped to independent values. We note here that the ability of a random
field to map nearby vectors to nearby values is restricted by the following
result [19]:

######  Theorem 2.

Let $f(t)$ be the characteristic function of a non-degenerate distribution.
Then there are positive constants $\alpha,\beta$ such that $f(t)\leq 1-\alpha
t^{2}$ for $|t|\leq\beta$.

Intuitively speaking, this result limits the capability of a random field to
“compress” the images of nearby vectors.

A natural candidate for a random field is the _Gaussian Random Field_(GRF), a
key construct in probability, applied mathematics, and physics [1, 5]. The
Gaussian random field (Definition 5) is also obtaining – up to a constant –
the best “compressibility” possible.

###  3.1. Gaussian Random Fields

######  Definition 5.

A (stationary) Gaussian Random Field on a set $D\subset\mathbb{R}^{d}$ is a
class of random real-valued functions on $D$, satisfying the following:
For every $x_{1}\ldots x_{n},x_{i}\in D$, the vector $(f(x_{1})...f(x_{n}))$
obeys a normal distribution with zero mean and a covariance matrix $C$ given
by $C_{i,j}=K(x_{i},x_{j})$, for some positive definite kernel $K(,)$.

Recall that a kernel is positive definite if it satisfies
$\int_{D}K(x,y)f(x)f(y)dxdy\geq 0$ for every $f\in L_{2}(D)$; this also means
that for every set of vectors $x_{1}\ldots x_{n}$ in $D$, the matrix
$K(x_{i},x_{j})$ is positive definite. Since our only interest is in
_stationary and isotropic_ kernels – i.e., which depend only on $\|x-y\|$ – we
make the standard assumption that $K(x,y)\triangleq K\left(\|x-y\|\right)$.

The simplest – yet very important – GRF, is the one defined by

$K(t)=\exp\left(-\lambda\|t\|^{2}\right)\text{ , }$

for some constant $\lambda>0$. In this case, for each $x\in D$, $f(x)$ is a
standard normal random variable, and the covariance matrix of
$(f(x_{1})...f(x_{n}))$ is given by
$C_{i,j}=\exp\left(-\lambda\|x_{i}-x_{j}\|^{2}\right)$. Note that there are
two Gaussian functions involved: the first defines the covariance matrix of
the second normal distribution. We shall refer to this field as the _basic
GRF_. In Sections 6 and 7, more general GRF’s will be studied.

As demonstrated in Figure 4, for both one and two dimensional fields, as
$\lambda$ increases, the field become more “rough”, and the correlation
between its values at $x,y$ drops faster as $\|x-y\|$ increases.

Figure 4. Examples of Gaussian random fields. Top: 1D field, $\lambda=100$.
Second: 1D field, $\lambda=10000$. Third: 2D field, $\lambda=100$. Bottom: 2D
field, $\lambda=10000$.

Once a GRF is defined, we can use its _maximal value_ on $S$ in order to
decide whether property $\Pi_{\text{YES}}$ or $\Pi_{\text{NO}}$ holds, as will
be described next. We note that if $S$ is bounded (e.g. in a ball of fixed
radius), the maximal value remains bounded even _as the cardinality of $S$
tends to infinity_, in stark contrast to the FM hashing, for which the
estimate of the cardinality (i.e. the inverse of the minimal value) will tend
to infinity.

The proposed solution searches for the optimal values of the GRF parameter
$\lambda$, as well as that of a threshold $T$, which maximize the following
for a random $g$ drawn from the basic GRF (where $\max(g)$ always refers to
the maximal value over $S$):

$\Pr\left(\max(g)\geq
T\mid\text{case}\;\Pi_{\text{NO}}\right)-\Pr\left(\max(g)\geq
T\mid\text{case}\;\Pi_{\text{YES}}\right)$ (1)

Noting that $\Pi_{\text{NO}}$ and $\Pi_{\text{YES}}$ are satisfied by a large
class of sets $S$, in this work we follow a “worst case scenario”, and
maximize the expression in Equation (1) under the following assumptions:

  * •

Case $\Pi_{\text{YES}}$: $S$ consists of $k_{1}$ balls of radius
$\varepsilon$, i.e. it is _maximal_ (with respect to set containment) among
all sets satisfying $\Pi_{\text{YES}}$. Further, we assume no lower bound on
the pairwise distances between these $k_{1}$ balls.

  * •

Case $\Pi_{\text{NO}}$: $S$ consists of exactly $k_{2}$ vectors, which are
equidistant from each other, all distance equal to $\delta$. Since we seek a
general bound on the field’s maximal value, we will ignore the fact that this
condition cannot always hold; for example, there are no 4 equidistant vectors
in $\mathbb{R}^{2}$.

Given the optimal parameter values, the decision whether $\Pi_{\text{YES}}$ or
$\Pi_{\text{NO}}$ holds proceeds by estimating Equation (1), by drawing a
fixed number of random fields, and computing the empirical probability of
their maxima to exceed $T$. Note that if $S$ is distributed among nodes, it
suffices for each node to send an estimate of the maximum of the GRF on its
part of $S$ to a coordinator node, which then simply takes the maximum of
these local maxima; thus the algorithm is very appropriate for a distributed
and/or dynamic setting, as the communication overhead does not depend on the
cardinality of $S$.

There remains the question of how to compute the two probabilities in Equation
(1), which we address next. One technical tool we use is the Slepian Lemma,
for GRFs.

######  Theorem 3.

The Slepian Lemma [24]: let
$(X_{1}\ldots X_{n})$ and $(Y_{1}\ldots Y_{n})$ be normal random variables
with zero expectation, and further satisfying $E[X_{k}^{2}]=E[Y_{k}^{2}]$ for
every $k$, and $E[Y_{k}Y_{\ell}]\geq E[X_{k}X_{\ell}]$ for every $k,\ell$.
Then, for every $T$, $\Pr(\max{X_{k}}\geq T)\geq\Pr(\max{Y_{k}}\geq T)$.

That is, if the expectation and variance of $X_{k},Y_{k}$ are equal for every
$k$, and all $X_{k},X_{\ell}$ pairs are “less correlated” than
$Y_{k},Y_{\ell}$, then the probability of $\max{X_{k}}$ to exceed $T$ is
larger than that that of $\max{Y_{k}}$ to exceed $T$, for every $T$.
Surprisingly, this intuitive result is _not_ true for all types of random
variables.

###  3.2. Probability for the Maximum of the GRF to Exceed $T$ on a set
Containing $k_{2}$ vectors with all Mutual Distances $\geq\delta$.

This problem reduces to that of computing the so-called exceedance probability
for a $k_{2}$-dimensional normal distribution with zero expectation and a
covariance matrix $C$ satisfying $C_{k,k}=1$ for all $k$, and
$k\neq\ell\rightarrow C_{k,\ell}\leq\exp\left(-\lambda\delta^{2}\right)$. The
following theorem allows to bound this probability from below by solving for
the case $C_{k,k}=1,C_{k\neq\ell}=\exp\left(-\lambda\delta^{2}\right)$:

Due to the Slepian Lemma, we need only to compute the probability for the
maximum to exceed $T$ for a normal vector $(X_{1}\ldots X_{k_{2}})$ with zero
expectation and covariance matrix
$C_{k,k}=1,C_{k\neq\ell}=\exp\left(-\lambda\delta^{2}\right)$. This
probability can be computed as follows: denote
$\exp\left(-\lambda\delta^{2}\right)$ by $\rho$. The random vector
$(X_{1}\ldots X_{k_{2}})$ can be modeled by

$x_{i}=\sqrt{\rho}N_{0}+\sqrt{1-\rho}N_{i},i=1\dots k_{2},$

where $N_{0},N_{1}\ldots N_{k_{2}}$ are i.i.d standard normal random variables
(it is easy to verify that each $X_{i}$ is a standard normal random variable,
and that $i\neq j\rightarrow\mathrm{Cov}\left(X_{i},X_{j}\right)=\rho$). For
any $T_{1}$, the probability density of
$\max{\\{\sqrt{1-\rho}N_{1},\ldots,\sqrt{1-\rho}N_{k_{2}}\\}}$ to equal
$T_{1}$ is given by

$\frac{\sqrt{2}\,{\mathrm{e}}^{-\frac{\mathit{T_{1}}^{2}}{2\left(1-\rho\right)}}}{2\sqrt{\pi}\,\sqrt{1-\rho}}\left(\frac{\mathrm{erf}\\!\left(\frac{\sqrt{2}\,\mathit{T_{1}}}{2\sqrt{1-\rho}}\right)}{2}\right)^{k_{2}-1},$

and the probability for the maximum of $x_{i}$ to exceed $T$ is therefore
given by

$\displaystyle\int_{-\infty}^{\infty}\left(\frac{\sqrt{2}\,{\mathrm{e}}^{-\frac{x^{2}}{2\rho}}}{2\sqrt{\pi}\,\sqrt{\rho}}\right)\cdot\,\frac{\sqrt{2}\,{\mathrm{e}}^{-\frac{\mathit{T-x}^{2}}{2\left(1-\rho\right)}}}{2\sqrt{\pi}\,\sqrt{1-\rho}}\left(\frac{\mathrm{erf}\\!\left(\frac{\sqrt{2}\,\mathit{T-x}}{2\sqrt{1-\rho}}\right)}{2}\right)^{k_{2}-1}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!dx,$

which can be evaluated numerically.

###  3.3. Probability for the Maximum of the GRF to Exceed $T$ on $k_{1}$
Balls with Radius $\varepsilon$ each

There are no known closed-form expressions for this probability, even in the
one-dimensional case (i.e. an interval), and the probability is computed using
a sampling approach, or approximations [4, 20]. In Section 4 we define a
different field, for which there is a closed-form expression for the
exceedance probability over a ball in any dimension. In order to derive an
absolute bound, we treat the field’s values on the $k_{1}$ balls as
independent, which produces the bound $1-\left(1-\Pr\left(\text{exceedance for
a single ball}\right)\right)^{k_{1}}$.

###  3.4. Examples for Optimal Parameters $T,\lambda$

Next follow some examples for choosing the optimal hyper-parameters
$T,\lambda$, including a case in which the proposed paradigm fails, i.e. no
$\lambda,T$ exist such that the exceedance probability for ”non-clusterable”
($\Pi_{\text{NO}}$) is higher than for ”clusterable” ($\Pi_{\text{YES}}$).

####  Example 1

Here, we seek to differentiate between case $\Pi_{\text{YES}}$ being one
cluster of size 0.01 (i.e. $k_{1}=1,\varepsilon=0.005$, and case
$\Pi_{\text{NO}}$ corresponding to a set which contains two vectors at least
0.02 apart ($k_{2}=2,\delta=0.02$), in dimension $d=1$. In Figure 6, the
values of Equation (1) are plotted for a range of $\lambda,T$. The optimal
parameters are obtained at the highest point on the surface.

Looking at $\lambda$ cross-sections of the 2D function in Figure 6 allows an
intuitive understanding of the dependence on $\lambda$. If $\lambda$ is very
small, then the covariance between $f(x),f(y)$, $\exp(-\lambda(x-y)^{2})$, is
very close to 1, hence the field is nearly constant, and its maximum is nearly
identical for cases $\Pi_{\text{YES}}$ and $\Pi_{\text{NO}}$. At the other
extreme, where $\lambda$ is very large, the field’s behavior approaches that
of white noise, and the maximum is nearly independent of the extent of $S$
(while it does depend on the number of vectors in $S$). For “good” $\lambda$
values, e.g. 3000, the correlation for case $\Pi_{\text{YES}}$ is bounded from
below by $\exp(-3000\cdot 0.01^{2})=0.741$, while for case $\Pi_{\text{NO}}$
it reaches $\exp(-3000\cdot 0.02^{2})=0.301$. This difference results in the
exceedance probability of the field in case $\Pi_{\text{NO}}$ being higher
than in case $\Pi_{\text{YES}}$, even when the number of vectors in case
$\Pi_{\text{YES}}$ is unbounded – as we compute the exceedance probability for
a field on an interval of length 0.01.

#### Example 2

There are however cases in which the proposed paradigm can fail, as there are
no values of $\lambda,T$ for which Eq. (1) obtains a positive value. This is
due to the fact that for case $\Pi_{\text{YES}}$, the exceedance probability
must be bounded on a continuous set, as opposed to case $\Pi_{\text{NO}}$,
where the field can only take values at a (typically small) finite set. Such
an example is when $\Pi_{\text{YES}}$ is defined by
$k_{1}=3,\varepsilon=0.05$, and $\Pi_{\text{NO}}$ by $k_{2}=4,\delta=0.1$.
Since the difference of the exceedence probabilities (Equation (1)) is always
negative (Figure 6), there are no parameters $\lambda,T$ which can be applied.

Figure 5. An example of optimal $\lambda,T$ parameters

Figure 6. An example in which the proposed paradigm fails

###  3.5. Constructing the Gaussian Random Field

There are no closed-form expressions for constructing a GRF, even in dimension
1. Here we propose an efficient approximation, which can be extended to higher
dimensions. The field is restricted to a finite interval, which without loss
of generality we assume to be the unit interval $[0,1]$.

Define

$h(x)=\displaystyle\sum_{k=0}^{n}N_{k}(0,1)w_{k}\cos(\pi k(x-t)),\;\;t\sim
U[-1,1]$ (2)

Since $h$ is the sum of normal variables, it is normal. Also, as can be
verified by direct integration, for any $x,y\in[0,1]$ the following holds:

$\mathrm{Cov}\left(h(x),h(y)\right)=\frac{1}{2}\displaystyle\sum_{k=0}^{n}w_{k}^{2}\cos\left(\pi
k(x-y)\right)$

Recall that for a GRF with parameter $\lambda$,
$\mathrm{Cov}(h(x),h(y))=\exp\left(-\lambda(x-y)^{2}\right)$. Therefore, we
seek a valid approximation on the interval $[0,1]$:

$\exp\left(-\lambda t^{2}\right)\approx\sum_{k=0}^{n}\beta_{k}\cos(\pi kt)$

Such an approximation can easily be found – it is a truncated Fourier
expansion of the Gaussian. However, there is a caveat – the Fourier
coefficients are not all positive, hence cannot be expressed as $w_{k}^{2}$.
While an optimal approximation using only positive coefficients can be sought,
for our purpose it is enough to simply omit the negative coefficients from the
expansion. Since over the continuous domain, the Fourier transform of a
Gaussian is also a Gaussian, all its coefficients are positive.

This idea can be extended to higher dimensions $d$; alas, the size of the
corresponding approximation increases with $d$ and $\lambda$ increases . The
reason for the dependence on $\lambda$ is that more Fourier coefficients are
required to approximate a Gaussian with a small standard deviation than with a
large one. The number of terms required to approximate a $d$-dimensional
Gaussian behaves as $\mathcal{O}(\lambda^{d/2})$.

Figure 7. The covariance of a 1D GRF with $\lambda=300$ (blue) and the
covariance function of an approximation with 15 Fourier elements (red).

Figure 8. The covariance of a 1D GRF with $\lambda=300$ (blue) and the
covariance function of an approximation with 25 Fourier elements (red).

###  3.6. Algorithm and Analysis

We now present our algorithm for distinguishing between the two cases,
$\Pi_{\text{YES}}$ and $\Pi_{\text{NO}}$. We assume that a random field was
chosen, based on the parameters $\varepsilon,\delta$ (see Sections 4, 6, and
7).
Input. Parameters: $\varepsilon,\delta\in\mathbb{R}_{+}$ and
$k_{1},k_{2}\in\mathbb{N}$. (1) Given the parameters
$\varepsilon,\delta,k_{1}$ and $k_{2}$ , find the $T$ and $\lambda$ 222This
step is done numerically. This is a main limitation of our method, and we
leave a more rigorous analysis of this step for future research. that maximize
Equation (1): $\Pr\left(\max(g)\geq
T\mid\text{case}\;\Pi_{\text{NO}}\right)-\Pr\left(\max(g)\geq
T\mid\text{case}\;\Pi_{\text{YES}}\right)\text{.}$ Denote: •
${\mathcal{C}}=\Pr\left(\max(g)\geq T\mid\text{case}\;\Pi_{\text{YES}}\right)$
• ${\mathcal{M}}=\Pr\left(\max(g)\geq
T\mid\text{case}\;\Pi_{\text{NO}}\right)$ (2) If no $\lambda,T$ exist for
which the equation above is strictly positive, return FAIL (3) Let
${\mathcal{P}}\leftarrow\Pr\left(\max(g)>T\right)$ over $S$333The probability
is computed empirically by generating a number of fields., for a set of
vectors $S\subseteq\mathbb{R}^{d}$. (4) If
${\mathcal{P}}<\frac{{\mathcal{M}}+{\mathcal{C}}}{2}$ return
$\Pi_{\text{YES}}$; otherwise, output $S\in\Pi_{\text{NO}}$.

Using standard probability techniques, the success rate of the algorithm can
be estimated.

###  3.7. Distributed Systems and Privacy

As for the FM algorithm, the approach described here for deciding between
$\Pi_{\text{YES}}$ and $\Pi_{\text{NO}}$ is very suitable in a distributed,
dynamic system: since the outcome depends only on the maximal value of the
field, it requires only a constant volume of communication per node, and is
also easy to adapt to the case of distributed data streams. The nodes merely
need to agree on a seed for generating the random field (as the same fields
must, of course, be applied on the entire data).

Interestingly, the above observation still allows a measure of privacy. Assume
for example a system with two nodes $N_{1},N_{2}$. If the maximal values on
the (identical) fields used by the nodes are well-correlated, $N_{2}$ can
infer on $N_{1}$’s data (for example, if $\Pi_{\text{YES}}$ holds, $N_{2}$
knows that $N_{1}$’s data can be covered by balls which also cover its own
data). However, if the maxima are uncorrelated, such inference is much more
difficult and partial, as the average values of the GRF are invariant to data
rotation and translation; therefore $N_{2}$ can only infer on $N_{1}$’s data
up to rotation and translation. We leave further analysis of this subject for
future research.

###  3.8. Dimensionality Reduction

The proposed algorithm is susceptible to curse of dimensionality, as the
maximal value of the GRF on a ball with a fixed radius increases with the
dimension, hence the exceedance probability for $\Pi_{\text{YES}}$ increases
with $d$, while for $\Pi_{\text{NO}}$ it does not. This is also true for the
other random fields we define (Sections 4, 7). As opposed to other problems,
ours cannot be directly remedied by applying random projections to low
dimensions (in the spirit of the famous Johnson-Lindenstrauss Lemma), since,
while preserving the ratio between distance, the absolute distances decrease;
however, the projection of a ball is a ball with the same radius. So while the
dimension decreases, so does $\delta$, making condition $\Pi_{\text{NO}}$ more
difficult to satisfy. The trade-off between reducing $d$ and reducing $\delta$
is an interesting topic for future research.

##  4\. Efficient Construction of Random Fields

While the GRF kernel is very suitable to our problem, it is computationally
expensive to construct, especially in high-dimensions, and its applications
are typically restricted to dimension 3 and below [5, 8]. Further, it is
untrivial to compute its exceedance probability on a ball; even for the 1D
case (interval), there are no closed-form expressions.

We start the analysis of other fields with the well-known method of Locally
Sensitive Hashing (LSH) [9]. The hash function, which accepts real valued
vectors as inputs, is defined (up to a scaling factor) by:

$h(x)=\left\lfloor{a\cdot x+b}\right\rfloor,\;a_{i}\sim N(0,\sigma),b\sim
U[0,1]$ (3)

where the normal distribution $N(0,\sigma)$ can be replaced by another _Stable
Distribution_.

In order to compare this hashing scheme to ours, we estimate the mean squared
distance between the hashed images of two vectors $x,y\in\mathbb{R}^{d}$. We
ignore the rounding operation (which is a necessity imposed by the requirement
to compare hash values). Direct integration yields:

$\displaystyle\mathop{\mathbb{missing}}{E}[(h(x)-h(y))^{2}]={(2\pi\sigma^{2})^{-\frac{d}{2}}}\cdot$
$\displaystyle\int_{-\infty}^{\infty}\ldots\int_{-\infty}^{\infty}\exp\left(-\frac{\|a\|^{2}}{2\sigma^{2}}\right)\left(a\cdot
x-a\cdot y\right)^{2}da_{1}\ldots da_{d}$
$\displaystyle=\sigma^{2}\|x-y\|^{2}$

It’s also easy to verify that $\mathop{\mathbb{missing}}{E}[f(x)]=0$, and
$E[f(x)f(y)]=x\cdot y$. Hence, the covariance of $h(x),h(y)$ does not tend to
0 when $\|x-y\|$ increases, as the GRF does; further, the maximal value of the
hash function will not be related to the number of clusters, but to the
diameter of the set. Hence, while LSH is highly successful in solving other
problems (e.g. nearest neighbor), it is not suitable to our goal.

###  4.1. Random Trigonometric Fields

In Section 3.5, we looked at a method to approximate a GRF using sums of
trigonometric functions, with random coefficients. Not surprisingly, the
larger $\lambda$ is, the more terms are required to approximate a field with
parameter $\lambda$, which satisfies
$\mathrm{Cov}\left(f(x),f(y)\right)=\exp(-\lambda\|x-y\|^{2})$, as the fields
become ”rougher” when $\lambda$ increases. It is natural to ask whether
similar fields can be obtained while using a smaller number of trigonometric
terms; it turns out that this can be achieved with one term, if the choice not
of the coefficients, but the argument of the trigonometric functions (i.e. the
frequency), is random 444In the _cosine field_ , [1], the frequency is
constant..

######  Definition 6 (Random Sine Field).

Let $x\in\mathbb{R}^{d}$. Define a random field – which will be referred to
hereafter as the _Random Sine Field_(RSF), with parameter $\sigma$, by:

$f(x)=\sin\left(a\cdot x+b\right)$ (4)

where $a\in\mathbb{R}^{d}$ is a random vector with components which are i.i.d
normal variables with zero expectation and $\sigma$ standard deviation, and
$b$ is uniformly distributed between $-\pi$ and $\pi$.

######  Lemma 4.

  1. (1)

for every $x\in\mathbb{R}^{d}$, $f(x)$ is a random variable with zero
expectation and a variance equal to $\frac{1}{2}$.

  2. (2)

The covariance of $f(x)$ and $f(y)$ is equal to:

$\frac{1}{2}\exp\left(-\frac{\sigma^{2}\|x-y\|^{2}}{2}\right).$

###### Proof.

We start with the expectation:

$\displaystyle
E[f(x)]=\frac{1}{2\pi}\int_{-\infty}^{\infty}\ldots\int_{-\infty}^{\infty}\frac{1}{(2\pi\sigma^{2})^{\frac{d}{2}}}$
$\displaystyle\int_{-\pi}^{\pi}\exp\left(-\frac{\|a\|^{2}}{2\sigma^{2}}\right)\sin(a\cdot
x+b)dbda_{1}\ldots da_{d}=0$

as follows immediately by first integrating over $b$.

For the variance, the proof also follows by direct integration:

$\displaystyle
E[f(x)f(y)]=\frac{1}{2\pi}\int_{-\infty}^{\infty}\ldots\int_{-\infty}^{\infty}\frac{1}{(2\pi\sigma^{2})^{\frac{d}{2}}}$
$\displaystyle\int_{-\pi}^{\pi}\exp\left(-\frac{\|a\|^{2}}{2\sigma^{2}}\right)\sin(a\cdot
x+b)\sin(a\cdot y+b)dbda_{1}\ldots da_{d}$
$\displaystyle=\frac{1}{2}\exp\left(\frac{-\sigma^{2}\|x-y\|^{2}}{2}\right)$

∎

The drawback of this field is that, as opposed to the basic GRF, it does not
satisfy the property that, for each $x_{1}\dots x_{n}$, the vector
$\left(f(x_{1})\ldots f(x_{n})\right)$ obeys a Gaussian distribution.
Therefore, we cannot a-priori use the Slepian Lemma [24], which guarantees
that, as the vectors in case $\Pi_{\text{NO}}$ move farther apart, the
exceedance probability of the maximum of the field increases. However, we were
able to prove an analogous result for the sine field:

######  Lemma 5.

The probability of the RSF to exceed any given value on a set of vectors, is
monotonically increasing in all the distances between the vectors 555We assume
that the vectors are normalized..

The proof of the Lemma is in section 5.

In addition to the very low computational complexity of generating the RSF –
even in high dimensions – it has another advantage: for a ball
$B_{\varepsilon}$ of radius $\varepsilon$ and a threshold $T$, there exists an
explicit, simple expression for the probability of the maximum of the field
over $B_{\varepsilon}$ to exceed $T$.

######  Lemma 6.

The probability of the RSF to exceed $T$ on $B_{\varepsilon}$ is given by

$\frac{\pi-2\arcsin(T)+2\varepsilon\sigma\sqrt{d}}{2\pi}$ (5)

where it is tacitly assumed that the numerator is $\leq 2\pi$.

###### Proof.

Clearly the probability of a random angle $\theta$ to satisfy
$\sin(\theta)\geq T$ is $\frac{\pi-2\arcsin(T)}{2\pi}$. It remains to compute
the probability that the image of $B_{\varepsilon}$ under a random hash
function intersects the interval of length $\frac{\pi-2\arcsin(T)}{2\pi}$ on
the unit circle in which $\sin(\theta)\geq T$. The image of $B_{\varepsilon}$
under $f(x)=a\cdot x$ is an interval of average length
$2\varepsilon\sigma\sqrt{d}$, as its endpoints are the inner products of $a$
with the vectors on $B_{\varepsilon}$ in the directions of $a$ and $-a$. The
average length of a $d$-vector whose coordinates are all drawn from
$N(0,\sigma)$ is $\sigma\sqrt{d}$. The probability for both intervals to
intersect is the sum of their lengths, hence the result follows. The geometric
idea behind the proof is provided in Fig. 9. ∎

Figure 9. Sketch for the proof of Lemma 6. The red ”dome” represents the range
of angles $\theta$ for which $\sin(\theta)\geq T$. The $d$-dimensional ball
with radius $\epsilon$ is projected onto a range of $\theta$ of length
$2\varepsilon\sigma\sqrt{d}$ (green). For the two intervals to intersect, the
starting point of the ”green” interval has to be between $\theta_{1}$ and
$\theta_{2}$.

Choosing optimal values for the parameters $T,\sigma$ proceeds as for the GRF
(Section 3), however they are easier to compute, since computing Eq. (1) is
simpler.

Fig. 10 provides a simple example, which captures the intuition behind the
ability of the RSF to distinguish conditions $\Pi_{\text{YES}}$ and
$\Pi_{\text{NO}}$: for small $\sigma$, the field is nearly constant, and no
differentiation between the aforementioned conditions is possible. When
$\sigma$ becomes too large, the field oscillates rapidly, and always obtains
the maximal value of 1 on a set of size $\varepsilon$, causing $Pr(\max(g)\geq
T\mid\Pi_{\text{YES}})$ to be larger than $Pr(\max(g)\geq
T\mid\Pi_{\text{NO}})$, as $\Pi_{\text{NO}}$ contains finitely many vectors.
The optimal value is obtained between these extreme values. This echoes the
choice of the optimal $\lambda$ for the GRF (Section 3); for small $\lambda$
the field is nearly constant, but as $\lambda$ increases, the field changes
more rapidly, until the expected maximum for $\Pi_{\text{YES}}$ will be larger
than for $\Pi_{\text{NO}}$. Intuitively, we seek values of $\lambda$ and
$\sigma$ for which the range of field values for $\Pi_{\text{NO}}$ is the
largest w.r.t to $\Pi_{\text{YES}}$.

Figure 10. An example of choosing the optimal $T,\sigma$ for the same
parameters as in Fig. 6 ($d=1,k_{1}=1,k_{2}=2,\delta=0.02,\varepsilon=0.01$).
As in Fig. 6, the difference between the exceedance probabilities of
conditions $\Pi_{\text{NO}}$ and $\Pi_{\text{YES}}$ is plotted, relative to a
range of values for $T,\sigma$. Note that in this case, RSF performs on par
with the GRF, as the maximal difference between $Pr(\Pi_{\text{NO}})$ and
$Pr(\Pi_{\text{YES}})$ is similar.

###  4.2. Random Fields and Stable Distributions

When defining the RSF, we randomly chose the vector $a$ such that each of its
components obeys a normal distribution. The choice of the distribution over
$a$ must be judicious, as demonstrated by the following example, which holds
even in the case of a 1D field: define a field by $f(x)=sin(ax+b)$, where
$b\sim U[-\pi,\pi]$ and $a\sim U[-A,A]$. It is easy to verify that
$\mathrm{Cov}(f(x),f(y))=\frac{\sin(A(x-y)}{2(x-y)}$ (i.e. a sinc function at
$x-y$). This is inappropriate for our goal, as the covariance is not positive
nor monotonically decreasing in $|x-y|$.

However, just as the normal distribution was used to construct the RSF in
Definition 4, we can use another well-known example of a stable distribution
for the coefficients of the vector $a$, which defined the field at $x$ by
$\sin(a\cdot x+b)$: $\Pr(a_{i})=\frac{k}{\pi(1+k^{2}a_{i}^{2})}$ (where $b$ is
uniformly distributed in $[-\pi,\pi]$, and $k>0$). Allowing a slight abuse of
terminology, we shall refer to it as the Stable Field. A somewhat more
elaborate computation (which will not be repeated here) than for the case of
the RSF, yields that in this case, the covariance matrix of the field is given
(up to a multiplicative constant) by:

$\mathrm{Cov}(f(x),f(y))=\exp\left(-\frac{\|x-y\|_{1}}{k}\right)$

while this field is also viable for our problem (it is stationary, and the
covariance monotonically decreases with the $L_{1}$ distance), its performance
will be inferior to that of the GRF and RSF; this is due to the fact that the
decrease in the covariance’s values from ”distance $\varepsilon$” to ”distance
$\delta$” will be slower. If e.g. $\delta=10\varepsilon$, and we set the
covariance to be 0.99 for distance $\varepsilon$, then for distance $\delta$
it will be 0.368 for the RSF and 0.905 for the stable field. Still, an
interesting topic for future work is the relation between stable distributions
and random fields, and their application to the problem studied here.

Figs. 12,12 include one and two dimensional examples of the stable field. Note
that these are not samples of the random field, which is just a sinusoidal
function, but a visualization of the covariance. These fields tend to be more
”jittery” than the GRF and RSF.

Figure 11. An example of a 1D stable field with $k=10$

Figure 12. An example of a 2D stable field with $k=10$

##  5\. Slepian-type result for RSF

In this section we prove an analogous version of the Slepian Lemma for the
sine random field. Recall that the “sine field” is defined by

$R(u)=\sin(a\cdot u+b),a_{i}\sim N(0,\sigma),b\sim U[-\pi,\pi]$

(assume hereafter $\sigma=1$, the proof remains the same for every $\sigma$).

As a direct computation proves 4, the covariance of $R(x),R(y)$ equals
$\exp(-||x-y||^{2})$. The decision whether properties $\Pi_{\text{NO}}$ or
$\Pi_{\text{YES}}$ hold for a set of points $S$ depends only on the exceedance
probability: given $T$, the probability that the maximum value $R()$ attains
on $S$ is at least $T$.

We next prove Lemma 5, which states that as the points of $S$ move farther
apart, the expectation of the maximum of the sine field increases.

Let us start with two points, $u_{1},u_{2}\in{\mathbb{R}}^{d}$. Assume
$||u_{i}||=1$. We then have, from standard results about Gaussian random
variables:

$R(u_{i})=\sin(x_{i}+b),x_{i}\sim N(0,1),\mathrm{Cov}(x_{1},x_{2})=u_{1}\cdot
u_{2}$

Since $u_{1}\cdot
u_{2}=\frac{1}{2}(||u_{1}||^{2}+||u_{2}||^{2}-||u_{1}-u_{2}||^{2})$ and
$||u_{i}||=1$, $\mathrm{Cov}(x_{1},x_{2})$ decreases when $||u_{1}-u_{2}||$
increases.

According to [25], when applied to
$f(x,y)=\int\limits_{-\pi}\limits^{\pi}\max(\sin(x+b),\sin(y+b))db$, it
suffices to prove that

$\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}p(x,y)\frac{\partial^{2}f}{\partial
x\partial y}dxdy\leq 0$ (6)

Next, we look at $\frac{\partial^{2}f}{\partial x\partial y}$. This is a
little tricky, as the mixed derivative of $\max(u,v)$ is 0 when $u\neq v$, and
does not exist when $u=v$. However, standard arguments allow to relate to it
as $-\frac{1}{2}\delta(u,v)$. Inserting the derivative into the integral which
defines $f(x,y)$, we have

$\frac{\partial^{2}f}{\partial x\partial
y}=\int\limits_{-\pi}\limits^{\pi}\cos(x+b)\cos(y+b)\frac{\partial^{2}\max(u,v)}{\partial
u\partial v}(u=\sin(x+b),v=\sin(y+b))db$

Therefore, the only points which contribute to the integrand are the values
$b$ for which $\sin(x+b)=\sin(y+b)$. When $x=y+2\pi k$, this holds for every
$b\in[-\pi,\pi]$, hence the integrand is a delta function at every point. In
addition, the integrand contains a $\cos(x+b)$ factor (and not
$\cos(x+b)\cos(y+b)$, as one of the cosines is canceled by the inner
derivative of the delta function). This means that along the lines $y=x+2\pi
k$, $f(x,y)$ is negative, and its value tends to infinity w.r.t relative to
the value at other points; therefore, the maximal value of the integral in Eq.
6 is obtained when the values of $p(x,y)$ obtain their smallest possible value
along the “ridge” $y=x$, i.e. when the covariance of $p(x,y)$ approaches -1,
and $p(x,y)$ approaches a degenerate Gaussian whose which equals 1 when
$y=-x$, and 0 otherwise.

For this Gaussian, the integral in Eq. 6 reduces to an integral over the line
$y=-x$. This one-dimensional integral can be readily computed by
distinguishing between the case $x=k\pi,y=-k\pi$ (in which the integrand is a
delta function at every point), and the other cases; its value turns out to be
0. Since this is the maximal value, the condition in Eq. 6 is fulfilled, and
the expectation of the maximum indeed increases when the covariance of $x,y$
decreases. The general case (more than two random variables) follows
immediately, since we need only look at all the partial derivatives involving
all pairs of variables.

##  6\. Lower Bounds and Optimal Kernels

We now look at a more general case, in which the covariance matrix of the
field is defined by any _positive definite kernel_ $k()$, that is, a function
from $\mathbb{R}^{d}$ to $\mathbb{R}$ with the property that, for every set of
vectors $x_{1}\ldots x_{n}$, the matrix
$C_{i,j}=k\left(\|x_{i}-x_{j}\|\right)$ is positive definite. This allows to
apply more general kernels, whose covariance matrix $C_{i,j}$ can be different
from the one of the basic GRF, for which $k(t)=\exp(-\lambda\|t\|^{2})$, and
$C_{i,j}=\exp(-\lambda\|x_{i}-x_{j}\|^{2})$. Many families of positive
definite kernels are known, for example the Matérn kernels [23], and
applications, for example in machine learning, can benefit from the choice of
a suitable kernel [12, 22]. Here we deal with the case $d=1$, and in Sections
7, 8 the general case is addressed.

As discussed in Section 3, an optimal kernel $k()$ for our problem is one
which, loosely speaking, treats vectors with distance smaller than
$\varepsilon$ as “nearly identical”, and vectors with distance larger than
$\delta$ as “nearly independent”. Ostensibly this can be achieved by applying
a “step function” kernel

$k(t)=\begin{cases}1&\text{ if }t\leq\alpha\\\ 0&\text{ if
}t>\alpha\end{cases}$

for some $\varepsilon<\alpha<\delta$. Alas, this kernel is not positive
definite, as can be seen even for a case with three vectors. Assume the
vectors are $x_{1},x_{2},x_{3}$, and satisfy
$\|x_{1}-x_{2}\|,\|x_{2}-x_{3}\|<\varepsilon$, and $\|x_{2}-x_{3}\|>\delta$.
Then the matrix defined by $C_{i,j}=k(\|x_{i}-x_{j}\|)$ is equal to

$\left(\begin{array}[]{ccc}1&1&0\\\ 1&1&1\\\ 0&1&1\end{array}\right)$

which is not positive definite, as its eigenvalues are
$1+\sqrt{2},1,1-\sqrt{2}<0$. Hence, it is impossible to realize it as the
covariance matrix of any field.

We next formalize the requirements for an optimal kernel function $k()$, given
the problem parameters $\varepsilon,\delta$. In order for the problem to be
well-defined, we normalize $k()$ by requiring $k(0)=1,k(\varepsilon)=c$. Since
$\varepsilon$ is the allowed size of a “tight” cluster, and we aim for the
field to map tight clusters to a very small range, we typically set $c$ to be
very close to 1, e.g. 0.99. We also require $k()$ to monotonically decrease,
since the correlation between the field values at $x,y$ must decrease as
$\|x-y\|$ increases. We also demand that $k()$ be non-negative, for the
following reason: since $k(0)=1$, if there exists a $t_{1}$ such that
$k(t_{1})<0$, then by continuity there must be a value $t_{2}$ such that
$0<t_{2}<t_{1}$, and $k(t_{2})=0$. But then, from the Slepian Lemma,

$\Pr\\{\max\\{f(0),f(t_{1})\\}\geq T\\}>\Pr\\{\max\\{f(0),f(t_{2})\\}\geq
T\\}$

in contradiction to the requirement that the exceedance probability of the
maximum increases as the vectors move farther apart.

All the requirements from $k()$ are now summarized:

  1. (1)

Domain: $k(t)$ is defined in the range $0\leq t\leq 1$.

  2. (2)

Normalization: $k(0)=1,k(\varepsilon)=c$.

  3. (3)

$k()$ is positive definite.

  4. (4)

$k()$ is positive.

  5. (5)

$k()$ is monotonically decreasing.

  6. (6)

Under conditions 1-5, $k(\delta)$ is minimal.

Kernels that satisfy 3,4 above are often referred to as _doubly positive_ ,
and are of considerable interest in physics and applied mathematics [14, 10,
13, 15]. From the celebrated Bochner Theorem [17], a kernel is doubly positive
iff both it and its Fourier transform are positive. The rate of decay of such
kernels was studied, unfortunately there do not exist concrete results for
solving condition 6.

A large family of candidate kernels is provided by the following theorem, due
to Polya [6]:

######  Theorem 7.

Let $\phi:\mathbb{R}\rightarrow[0,\infty]$ be even, decreasing, and convex.
Then it defines a positive definite kernel.

Alas, the convexity constraint severely hampers the goal of obtaining a
minimal value of $k(\delta)$. To illustrate this, observe that for
$\varepsilon=0.01,\delta=0.1,c=0.99$, the value of a convex kernel at 0.1 is
bounded from below by 0.9, while the Gaussian kernel (known to be doubly
positive), $\exp(-\lambda t^{2})$, with $\lambda$ chosen to satisfy
$k(0.01)=0.99$, satisfies $k(0.1)=0.366$.

The question therefore is: can the value $k(\delta)$ obtained by the basic GRF
(Gaussian kernel) be decreased further? To the best of our knowledge, there
are no known bounds on the minimal value. We now describe two methods for
bounding $k(\delta)$. Such bounds are very important, as they allow us to
determine how close a certain kernel is to the optimal one.

###  6.1. Bounding $k(\delta)$ using positive definiteness of kernel matrices

Let us define, for every non-negative integer $l$, $p_{l}=k(l\varepsilon)$. If
$k()$ is a positive definite kernel, then for every $x_{1},x_{2}\ldots x_{n}$,
the $n\times n$ matrix $M_{i,j}=k(\|x_{i}-x_{j}\|)$ is positive definite. This
allows to derive many inequalities which are useful in bounding $k_{l}$. The
most basic one is provided in the following simple lemma.

######  Lemma 8.

$p_{2}\geq 2p_{1}^{2}-1$.

###### Proof.

taking $x_{0}=1,x_{1}=\varepsilon,x_{2}=2\varepsilon$, the matrix

$M=\left[{\begin{array}[]{ccc}1&p_{1}&p_{2}\\\ p_{1}&1&p_{1}\\\
p_{2}&p_{1}&1\\\ \end{array}}\right]$

must be positive definite, hence its determinant is positive. But
$|M|=(1-p_{2})(1+p_{2}-2p_{1}^{2})$. Since $1-p_{2}\geq 0$, the result
follows. ∎

While useful, Lemma 8 does not always suffice to provide lower bounds on the
values of $p_{l}$ – actually, not even $p_{2}$. In order to see this, assume
$p_{1}=0.7$. Then, Lemma 8 provides the bound $p_{2}\geq-0.02$. Since $p_{2}$
must be positive, this provides a lower bound of 0; and indeed, the
corresponding kernel matrix

$\left[{\begin{array}[]{ccc}1&0.7&0\\\ 0.7&1&0.7\\\ 0&0.7&1\\\
\end{array}}\right]$

is positive definite. Now, the positivity and monotonicity constraints imply,
obviously, that $p_{3}=0$ as well. However, the corresponding kernel matrix,
constructed using the points $0,\varepsilon,2\varepsilon,3\varepsilon$ equals

$\left[{\begin{array}[]{cccc}1&0.7&0&0\\\ 0.7&1&0.7&0\\\ 0&0.7&1&0.7\\\
0&0&0.7&1\\\ \end{array}}\right]$

which is _not_ positive definite, having a strictly negative eigenvalue
($\approx-0.1326$).

Hence, the lower bound provided by Lemma 8 is non-optimal (unless
$\varepsilon>\frac{1}{3}$, in which case $3\varepsilon>1$ and is not a valid
point).

This non-optimality of this bound is due to its “restricted view”, so to say;
it is derived using only the points $0,\varepsilon,2\varepsilon$ and their
corresponding kernel matrices, while ignoring points with larger coordinates
(e.g. $3\varepsilon$).

The lesson of this simple analysis is that better lower bounds on $p_{l}$ can
be obtained by looking at larger kernel matrices, involving more and more
points. For example, let us look at a typical range of parameters:
$\varepsilon=0.01,c=0.99$, and try to bound $k(0.1)$ from below. In principle,
a bound can be found by solving the following optimization problem:

$\text{argmin }{p_{10}}\text{ subject to}\hskip 5.69054ptA_{i}\succeq
0,i=1,2...\ldots n.$

where $\succeq 0$ stands for positive definiteness, and $A_{i}$ is the kernel
matrix constructed using the values $p_{l}$ for $l=1\ldots i$. Alas,that will
result in very large kernel matrices, which renders the optimization problem
non-practical (for example, optimizing over determinants of matrices which are
larger than $10\times 10$ is very difficult, due both to the high degree of
the resulting polynomials and their sheer size). Moreover, the large matrices
encountered in the optimization process turn out to be extremely ill-
conditioned.

A more tractable solution is to use a hierarchy of bounds. For example,
$p_{10}$ can be bounded as follows: first, take $\varepsilon=0.01$, and bound
$p_{5}$. Then look at $\varepsilon=0.05$, and bound $p_{2}$ (which is
equivalent to bounding $p_{10}$ for $\varepsilon=0.01$).

###  6.2. Bounding $p_{5}$

In order to bound $p_{5}$, we looked at the following combinations of points:
$\\{0,\varepsilon,2\varepsilon\\}$,$\\{0,\varepsilon,2\varepsilon,3\varepsilon\\}$,
$\\{0,4\varepsilon,5\varepsilon\\}$, $\\{0,2\varepsilon,4\varepsilon\\}$. The
determinants of the resulting kernel matrices are:

$\displaystyle\left({p_{2}}-1\right)\left(2\,{{p_{1}}}^{2}-{p_{2}}-1\right),$
$\displaystyle\left({{p_{1}}}^{2}+2\,{p_{1}}\,{p_{2}}-{p_{1}}\,{p_{3}}+{{p_{2}}}^{2}-{p_{1}}-{p_{3}}-1\right)$
$\displaystyle\left({{p_{1}}}^{2}-2\,{p_{1}}\,{p_{2}}-{p_{1}}\,{p_{3}}+{{p_{2}}}^{2}+{p_{1}}+{p_{3}}-1\right),$
$\displaystyle
2\,{p_{1}}\,{p_{4}}\,{p_{5}}-{{p_{1}}}^{2}-{{p_{4}}}^{2}-{{p_{5}}}^{2}+1,$
$\displaystyle\left({p_{4}}-1\right)\left(2\,{{p_{2}}}^{2}-{p_{4}}-1\right)$

The relative simplicity of these expressions (polynomials of degree at most
three or products of thereof), allowed a closed-form solution using the
MapleTM software package, providing a lower bound of $0.756$ for $p_{5}$.

It remains to bound $p_{2}$, assuming $p_{1}=0.756$. The naive bound
$2p_{1}^{2}-1$ yields $0.15477$, however, when the points
$\\{3\varepsilon,4\varepsilon\\}$ and the corresponding kernel matrices are
also considered (that is, we bound $k(2\varepsilon)$ using also the values
$k(3\varepsilon),k(4\varepsilon)$), then a much better bound is obtained,
$0.284$.

We note here that, to the best of our knowledge, no better bounds exist;
however, ostensibly such bounds can be found, by refining the optimization
process.

We next provide a different approach for bounding the kernel values, using a
discrete approximation to $k()$. This approach typically produces better
bounds, and also allows to construct approximately optimal kernels.

###  6.3. Bounding the kernel values using a discrete Fourier approximation

In this approach, we derive a lower bound on $p_{l}$ by directly optimizing
over the kernel $k(t)$, while imposing requirements 1-6 described in this
section’s first part. Alas, optimization techniques such as the calculus of
variations lead to unsolvable equations (for example, the monotonicity
constraint is equivalent to imposing overall negativity on the derivative, a
problem which is generally intractable). Still, we can obtain good
approximations by using the common approach of discretizing $k(t)$, a
technique which is used in numerous problems. Here, we discretize both the
Fourier transform $K(u)$ and the spatial coordinate $t$, as follows.

  1. (1)

Define a Fourier approximation to $k(t)$ as
$\displaystyle\sum_{l=0}^{n}a_{l}\cos(l\pi t)$.

  2. (2)

Impose the condition $k(0)=1$ by the linear equality
$\displaystyle\sum_{l=0}^{n}a_{l}=1$.

  3. (3)

Impose the condition $k(\varepsilon)=c$ by the linear equality
$\displaystyle\sum_{l=0}^{n}a_{l}\cos(l\pi\varepsilon)=c$.

  4. (4)

Impose (approximate) monotonicity by demanding that the derivatives at a dense
grid, e.g. $\\{0,0.001\ldots 0.999\\}$, are all negative. Note that these are
linear inequalities in the sought coefficients $a_{l},l=0\dots n$. Exact
monotonicity can be directly checked on the resulting solution.

  5. (5)

Impose positive definiteness of the kernel by the conditions $a_{l}\geq
0,l=0\dots n$.

  6. (6)

Impose positivity of the kernel by $k(1)\geq 0$, another linear constraint
(this suffices, since $k(t)$ is monotonically decreasing).

  7. (7)

under the constraints in 2-6 above, minimize $k(\delta)$, which is, also, a
linear function of $a_{l}\geq 0,l=0\dots n$.

Thus, the problem of finding an optimal kernel for our purpose has been
approximated by a linear programming problem. Increasing $n$ and the grid over
which the condition $k^{\prime}(t)\leq 0$ is imposed, improves the
approximation.

Many experiments have indicated that in the ”normal” range of parameters, the
optimal kernel is nearly equal to the simpler Gaussian kernel; this is next
demonstrated.

###  6.4. Examples of optimal kernels

Below are two figures, comparing the basic GRF to two optimal kernels, for
given values of $\varepsilon,\delta,c$.

#### Fig. 14

The optimal and Gaussian kernels, for $\varepsilon=0.01,c=0.99,\delta=0.1$.
Note that the Gaussian kernel (blue) is very close to the optimal one (red),
although (as can be expected) the optimal kernel achieves a slightly lower
value at $\delta$.

#### Fig. 14

An example in which the optimal kernel is markedly better than the Gaussian:
$\varepsilon=0.3,c=0.6,\delta=0.7$. Its ”odd” shape notwithstanding, the red
curve depicts a viable, doubly positive kernel. Interestingly, note that sharp
”drop” in the kernel just before $\delta$, in which the optimal kernel
achieves a much lower value than the Gaussian one.

However, the choice of parameters in Fig. 14 is less appropriate for the
practical decision problems we aim to solve in this paper; and typical values
yield optimal kernels which are rather close to the basic Gaussian kernel.

Figure 13. Comparison of the optimal (red) vs. Gaussian kernels (blue), for
$\varepsilon=0.01,c=0.99,\lambda=0.1$.

Figure 14. An example of an optimal kernel (red) which substantially differs
from the Gaussian kernel (blue), for $\varepsilon=0.3,c=0.6,\lambda=0.7$.

##  7\. Optimal Kernels in Higher Dimensions

In this section, we extend the approach of Section 6.3 to more than one
variable, with the Fourier expansion being rvarepsilonlaced by its
multivariate generalization, the Dini series.

###  7.1. Statement of the Problem

Given a dimension $d\geq 1$ and parameters $c,\varepsilon\in[0,1]$, we denote
by $\mathcal{K}_{d}(c,\varepsilon)$ the class of kernels
$k:\mathbb{R}_{+}\to\mathbb{R}$ such that

  1. (1)

$k(t)\geq 0$ for all $t\geq 0$;

  2. (2)

$k$ is continuous and monotonically decreasing in the interval $[0,1]$;

  3. (3)

The map $x\in\mathbb{R}^{d}\mapsto k(\|x\|)$ is positive semidefinite666We
will often say that $k$ is positive semidefinite in $\mathbb{R}^{d}$., that
is, for any collection of vectors $\\{x_{j}\\}_{j=1}^{N}\subset\mathbb{R}^{d}$
the matrix $[k(\|x_{i}-x_{j}\||)]_{i,j=1,...,N}$ is positive semidefinite;

  4. (4)

$k(\varepsilon)=c$.

Then, given $\delta\in(\varepsilon,1)$, we want to minimize the quantity
$k(\delta)$ over $\mathcal{K}_{d}(c,\varepsilon)$. We let

$\kappa_{d}(c,\varepsilon,\delta)=\inf_{k\in\mathcal{K}_{d}(c,\varepsilon)}k(\delta).$

######  Definition 7.

We say that $k_{0}\in\mathcal{K}_{d}(c,\varepsilon)$ is a global minima for
$\delta$ if $k(\delta)\geq k_{0}(\delta)$ for all
$k\in\mathcal{K}_{d}(c,\varepsilon)$. In addition, we say $k$ is a unique if
$k(\delta)=k_{0}(\delta)$ implies $k=k_{0}$.

Classical results of Bochner [17] and Schoenberg [21] show that:

  1. (1)

These classes are nested:
$\mathcal{K}_{d+1}(c,\varepsilon)\subset\mathcal{K}_{d}(c,\varepsilon)$.

  2. (2)

$k:\mathbb{R}_{+}\to\mathbb{R}$ is positive semidefinite and $k(0)=1$ if and
only if

$k(\|x\|)=\int_{\mathbb{R}^{d}}e^{ix\cdot y}d\mu(y)$

for all $x\in\mathbb{R}^{d}$, where $\mu$ is some probability measure in
$\mathbb{R}^{d}$.

  3. (3)

If we let $\mathcal{K}_{\infty}(c,\varepsilon):=\bigcap_{d\geq
1}\mathcal{K}_{d}(c,\varepsilon)$, then $k$ must be of the form

$\displaystyle k(\|x\|)=\int_{0}^{\infty}e^{-\lambda\|x\|^{2}}d\mu(\lambda)$
(7)

for some probability measure $\mu$ on $\mathbb{R}_{+}$.

######  Proposition 9.

Let
$\kappa_{\infty}(c,\varepsilon,\delta)=\inf_{k\in\mathcal{K}_{\infty}(c,\varepsilon)}k(\delta)$.
Then

$\kappa_{\infty}(c,\varepsilon,\delta)=c^{(\delta/\varepsilon)^{2}}$

and the kernel $k(t)=e^{-t^{2}\varepsilon^{-2}\log(1/c)}$ is the unique global
minima.

###### Proof.

By Schoenberg’s results (7) we have that

$k(t)=\int_{0}^{\infty}e^{-\lambda t^{2}}d\mu(\lambda)$

for some probability measure $\mu$. In particular, since
$\delta/\varepsilon>1$, Jensen’s inequality gives

$c^{(\delta/\varepsilon)^{2}}=k(\varepsilon)^{(\delta/\varepsilon)^{2}}=\left(\int_{0}^{\infty}e^{-\lambda\varepsilon^{2}}d\mu(\lambda)\right)^{(\delta/\varepsilon)^{2}}\leq\int_{0}^{\infty}e^{-\lambda\delta^{2}}d\mu(\lambda)=k(\delta).$

Equality can only be attained above iff the function $\lambda\mapsto
e^{-\lambda\varepsilon^{2}}$ is constant in the support of $\mu$, that is,
$\mu$ is a Dirac delta at some point $\lambda_{0}$. Since $k(\varepsilon)=c$
we must have $\lambda_{0}=\varepsilon^{-2}\log(1/c)$. This completes the
proposition. ∎

###  7.2. Dini Series and the Parameter space for
$\mathcal{K}_{d}(c,\varepsilon)$

It will be convenient to assume that our data is contained in a ball of radius
$1/2$, hence we can restrict our attention to the values of $k$ for
$t\in[0,1]$ only. Moreover, in working with the Dini series the way it is set
up, once $k$ is determined via a _finite_ Dini series with nonnegative
coefficients in the interval $[0,1]$, we can easily extend $k$ to
$\mathbb{R}_{+}$ using the series itself, since each function
$x\in\mathbb{R}^{d}\mapsto A_{\nu}(j_{\nu+1,n}\|x\|)$ is positive definite.

The Dini series777Sometimes also called Fourier-Bessel series. of a function
$f:[0,1]\to\mathbb{R}$ is given by

$f(x)=\mathfrak{B}_{0}(x)+\sum_{n=1}^{\infty}c_{n}J_{\nu}(\lambda_{n}x),$ (8)

where $0<\lambda_{1}<\lambda_{2}<\ldots$ denote the positive zeros of the
function

$zJ_{\nu}^{\prime}(z)+HJ_{\nu}(z)=(H+\nu)J_{\nu}(z)-zJ_{\nu+1}(z).$ (9)

Here, $J_{\nu}$ is the Bessel function of the first kind of order
$\nu\geq-\frac{1}{2}$, and $H\in\mathbb{R}$ is some parameter. This Bessel
function arises naturally as the Fourier transform of the unit sphere in
$\mathbb{R}^{d}$. Indeed, letting

$A_{\nu}(z):=\Gamma(\nu+1)(z/2)^{-\nu}J_{\nu}(z),$

we have

$\displaystyle
A_{d/2-1}(\|x\|)=\frac{1}{\text{Area}(\mathbb{S}^{d-1})}\int_{\mathbb{S}^{d-1}}e^{ix\cdot
y}dy$ (10)

for all $x\in\mathbb{R}^{d}$. The Dini series which is most convenient for our
purposes is the one with $H=-\nu$, so that equation (9) is simply
$zJ_{\nu+1}(z)=0$ and $\\{\lambda_{n}\\}_{n\geq 1}$ are positive Bessel zeros
of $J_{\nu+1}(z)$. Applying the identity [27, 5.11-(8)]

$\int_{0}^{1}J_{\nu}(\alpha t)J_{\nu}(\beta t)tdt=\frac{\alpha
J_{\nu+1}(\alpha)J_{\nu}(\beta)-\beta
J_{\nu}(\alpha)J_{\nu+1}(\beta)}{\alpha^{2}-\beta^{2}}.$ (11)

for $\alpha=\lambda_{m}$ and $\beta=\lambda_{n}$ are distinct zeros of (9),
then one sees that this integral vanishes whenever $m\neq n$. One can then
rescale the results from [27, 18.33], and invoke the identity [27, 5.11-(11)],
$\int_{0}^{1}A_{\nu}^{2}(\lambda_{n}t){t^{2\nu+1}}dt=\frac{A_{\nu}^{2}(\lambda_{n})}{2},$
in order to derive the following proposition.

######  Proposition 10.

Let $\nu\geq-\frac{1}{2}$. For every measurable $f:[0,1]\to\mathbb{C}$ with

$\|f\|_{\nu}:=\int_{0}^{1}|f(t)|^{2}{t^{2\nu+1}}dt<\infty$

we have that

$f(x)=a_{0}+\sum_{n=1}^{\infty}a_{n}{A_{\nu}(\lambda_{n}t)}$ (12)

with convergence in the norm $\|f\|_{\nu}$, where $\\{\lambda_{n}\\}_{n\geq
1}$ denote the positive zeros of the Bessel function $J_{\nu+1}$,

$a_{n}=\frac{2}{A_{\nu}(\lambda_{n})^{2}}\int_{0}^{1}f(t)A_{\nu}(\lambda_{n}t){t^{2\nu+1}dt},$
(13)

for all $n\geq 1$, and
$a_{0}=\tfrac{1}{2(\nu+1)}\int_{0}^{1}f(t){t^{2\nu+1}dt}$. Moreover, if $f$ is
continuous and of bounded variation in $[0,1]$, then the Dini series (12) of
$f$ converges absolutely and uniformly in $[0,1]$.

######  Remark 1.

For now on we set $\nu=d/2-1$, $\\{j_{\nu,n}\\}_{n\geq
1}=\text{Zeros}(J_{\nu})\cap\mathbb{R}_{+}$ and
$f_{d,n}(t)=A_{\nu}(j_{\nu+1,n}t)$ for $n\geq 1$ (note that
$f_{1,n}(t)=\cos(\pi nt)$). Identity (10) shows that any of the functions
$x\in\mathbb{R}^{d}\mapsto f_{d,n}(\|x\|)$ is indeed positive definite.

It is then straightforward to deduce the following result.

######  Proposition 11.

Let $k:\mathbb{R}_{+}\to\mathbb{R}$ be a continuous, of bounded variation and
positive semidefinite in $\mathbb{R}^{d}$. Then its Dini series

$k(t)=a_{0}+\sum_{n\geq 1}a_{n}f_{d,n}(t),$

for $t\in[0,1]$, has only nonnegative coefficients, that is, $a_{n}\geq 0$ for
all $n\geq 0$.

Observing that for all $0<t<1$ we have $|f_{d,n}(t)|<f_{d,n}(0)=1$, we
conclude that if $k\in\mathcal{K}_{d}(c,\varepsilon)$ then

$\displaystyle k(t)=1-\sum_{n\geq 1}a_{n}(1-f_{d,n}(t)).$ (14)

Moreover, letting $g_{d,n}(t)=\frac{1-f_{d,n}(t)}{1-f_{d,1}(t)}$ we conclude
that

$\displaystyle
k(t)=1-(1-c)\frac{1-f_{d,1}(t)}{1-f_{d,1}(\varepsilon)}+{(1-f_{d,1}(t))}\sum_{n\geq
2}a_{n}(g_{d,n}(\varepsilon)-g_{d,n}(t)).$ (15)

An important thing to notice is that the function

$\tilde{k}(t)=1-(1-c)\frac{1-f_{d,1}(t)}{1-f_{d,1}(\varepsilon)}$

belongs to $\mathcal{K}_{d}(c,\varepsilon)$ iff

$c\geq\frac{f_{d,1}(\varepsilon)-f_{d,1}(1)}{1-f_{d,1}(1)}.$

To see this, it is easy to show using basic facts about Bessel functions that
$f_{d,1}(t)$ is decreasing and that $f_{d,1}(1)<0$, and so $k_{1}$ is
decreasing and the condition above is just $\tilde{k}(1)\geq 0$.

##  8\. Numerical Evaluation

###  8.1. Heuristic

Identity (15) shows that if $g_{d,n}(\varepsilon)>g_{d,n}(\delta)$ for all
$n\geq 2$ then

$k(\delta)\geq\tilde{k}(\delta),$

and the kernel $\tilde{k}(t)$ is the unique global optimal. Numerically, one
can see that the functions $g_{d,n}(t)$ decay fast when $n\geq 2$ and $n=O(d)$
(especially for large $d$). These functions are essentially decreasing for all
$t\in[0,1]$ and this indicates there is no improvement in taking $a_{n}>0$ for
$n=O(d)$. Moreover, for $n\gg d$, the functions $g_{d,n}$ start to oscillate
rapidly and it becomes very difficult to enforce the condition
$k^{\prime}(t)\leq 0$, so there is no advantage in setting $a_{n}>0$ for $n\gg
d$ either. This explains what we observe in our numerical experiments of
Section 8: for small $\varepsilon$ and large $\delta$, it is almost always the
case that the kernel $\tilde{k}(t)$ is optimal, since the monotonicity of the
$g_{d,n}$ sets $a_{n}\approx 0$ for $n=2,..,O(d)$ and the condition
$k^{\prime}(t)\leq 0$ forces the remaining $a_{n}$ to vanish.

Let $\mathop{\mathbb{missing}}{E}_{d}(c)$ be the inverse of the function
$\varepsilon\mapsto c=\frac{f_{d,1}(\varepsilon)-f_{d,1}(1)}{1-f_{d,1}(1)}$.
Note that $c\mapsto\mathop{\mathbb{missing}}{E}_{d}(c)$ is decreasing with
$\mathop{\mathbb{missing}}{E}_{d}(0)=1$ and
$\mathop{\mathbb{missing}}{E}_{d}(1)=0$. In particular, given $c$, the kernel

$\tilde{k}(t)=1-(1-c)\frac{1-f_{d,1}(t)}{1-f_{d,1}(\varepsilon)}$

belongs to $\mathcal{K}_{d}(c,\varepsilon,\delta)$ exactly when
$\varepsilon\geq\mathop{\mathbb{missing}}{E}_{d}(c)$. For most practical
purposes we set $c\approx 0.9$. In the following table we estimate some values
of $\mathop{\mathbb{missing}}{E}_{d}(c)$ for $c=0.9$

$d$ | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10
---|---|---|---|---|---|---|---|---|---|---
$\mathop{\mathbb{missing}}{E}_{d}(0.9)$ | 0.2048 | 0.1991 | 0.1938 | 0.1890 | 0.1845 | 0.1803 | 0.1764 | 0.1728 | 0.1694 | 0.1661

It appears that the function $\mathop{\mathbb{missing}}{E}_{d}(c)$ decreases
with $d$ for every fixed $c$. This table shows that in low dimensions taking
$\varepsilon=0.21$ for a $0.9\leq c\leq 1$ guarantees that
$\tilde{k}(t)\in\mathcal{K}_{d}(c,\varepsilon,\delta)$. We have performed
numerous numerical experiments for $d=1,...,10$,
$\varepsilon\gtrapprox\mathop{\mathbb{missing}}{E}_{d}(c)$, several different
values of $c\lessapprox 1$ and $\delta>\varepsilon$, using the following
setup.

Linear Programming Setup. Given $N,M\geq 1$
$\displaystyle\begin{split}\text{Minimize}\quad\sum_{n=0}^{N}a_{n}f_{d,n}(\delta)\quad:\quad&a_{n}\geq
0\ \ \text{for each }n=0,...,N;\quad\sum_{n=0}^{N}a_{n}=1;\\\
&\sum_{n=0}^{N}a_{n}f_{d,n}(\varepsilon)=c;\quad\sum_{n=0}^{N}a_{n}f_{d,n}(1)\geq
0;\\\ &\sum_{n=1}^{N}a_{n}f_{d,n}^{\prime}(m/M)\leq 0\ \ \text{for each
}m=0,...,M.\end{split}$ (16)

Derivatives can be computed by the formula
$f_{d,n}^{\prime}(t)=-j_{\nu+1,n}^{2}tA_{\nu+1}(j_{\nu+1,n}t)/d$. Good results
only appear when we take $M>N^{2}$. In all cases, the kernel $\tilde{k}$ seem
to be the optimal choice, confirming our heuristic assumption.

######  Conjecture 12 (Heuristics).

If $c\lessapprox 1$ and
$\varepsilon\gtrapprox\mathop{\mathbb{missing}}{E}_{d}(c)$ then the kernel

$\tilde{k}(t)=1-(1-c)\frac{1-f_{d,1}(t)}{1-f_{d,1}(\varepsilon)}$

is the unique global minima of $\mathcal{K}_{d}(c,\varepsilon)$ for any
$\delta\in(\varepsilon,1]$.

The extreme case when $\varepsilon\ll\mathop{\mathbb{missing}}{E}_{d}(c)$
seems to be quite challenging, even in the desired range of parameters such as
$c\approx.9$ and $\delta\approx 1$. It seems that one can actually make
$k(t)=0$ for $t\geq\delta$ at the expense of $\varepsilon$ being very small,
see Figure 15.

$0$$0.1$$0.3$$0.5$$0.7$$0.9$$1$$0.2$$0.4$$0.6$$0.8$$1$$(d,c,\varepsilon,\delta)=(1,0.99,0.01,0.10)$$0$$0.1$$0.3$$0.5$$0.7$$0.9$$1$$0.2$$0.4$$0.6$$0.8$$1$$(d,c,\varepsilon,\delta)=(1,0.60,0.30,0.70)$$0$$0.1$$0.3$$0.5$$0.7$$0.9$$1$$0.2$$0.4$$0.6$$0.8$$1$$(d,c,\varepsilon,\delta)=(5,0.90,0.10,0.70)$$0$$0.1$$0.3$$0.5$$0.7$$0.9$$1$$0.2$$0.4$$0.6$$0.8$$1$$(d,c,\varepsilon,\delta)=(5,0.99,0.001,0.10)$$0$$0.1$$0.3$$0.5$$0.7$$0.9$$1$$0.2$$0.4$$0.6$$0.8$$1$$(d,c,\varepsilon,\delta)=(10,0.90,0.10,0.70)$$0$$0.1$$0.3$$0.5$$0.7$$0.9$$1$$0.2$$0.4$$0.6$$0.8$$1$$(d,c,\varepsilon,\delta)=(10,0.99,0.01,0.10)$
Figure 15. Examples of kernels for varying dimensions and parameters.

##  9\. Conclusions

We have presented a probabilistic hashing method, based on random fields,
which allows, for some parameter values, to differentiate between
”clusterable” and ”non-clusterable” data. The method easily extends to
distributed, dynamic environments, requiring only a constant communication
overhead, and maintaining a degree of privacy.

Future work will study further field types, as well as solving decision
problems for other geometric properties (e.g intrinsic dimensionality). We
also plan to study ”softer” definitions of properties $M,C$ (e.g. which still
hold if a certain percentage of the data is modified).

## References

  * [1] R. J. Adler and J. E. Taylor.  Random fields and geometry.  Springer Science & Business Media, 2009.
  * [2] N. Alon, S. Dar, M. Parnas, and D. Ron.  Testing of clustering.  SIAM Rev., 46(2):285–308, 2004.
  * [3] N. Alon, Y. Matias, and M. Szegedy.  The space complexity of approximating the frequency moments.  In G. L. Miller, editor, Proceedings of the Twenty-Eighth Annual ACM Symposium on the Theory of Computing, Philadelphia, Pennsylvania, USA, May 22-24, 1996, pages 20–29. ACM, 1996.
  * [4] J.-M. Azaïs and M. Wschebor.  A general expression for the distribution of the maximum of a gaussian field and the approximation of the tail.  Stochastic Processes and their Applications, 118(7):1190–1218, 2008.
  * [5] R. Bansal, D. X. Lawrence H. Staib, H. Zhu, and B. S. Peterson.  Statistical analyses of brain surfaces using gaussian random fields on 2-d manifolds.  IEEE Transactions on Medical Imaging, 26:46–57, 2006.
  * [6] C. Berg and G. Forst.  Potential theory on locally compact abelian groups.  In Potential Theory on Locally Compact Abelian Groups, pages 8–30. Springer-Verlag Berlin Heidelberg, 1975.
  * [7] S. Chakraborty, R. Pratap, S. Roy, and S. Saraf.  Helly-type theorems in property testing.  In A. Pardo and A. Viola, editors, LATIN 2014: Theoretical Informatics, pages 306–317, Berlin, Heidelberg, 2014. Springer Berlin Heidelberg.
  * [8] P. E. Creasey and A. Lang.  Fast generation of isotropic gaussian random fields on the sphere.  Monte Carlo Methods and Applications, 24(1):1–11, 2018.
  * [9] M. Datar, N. Immorlica, P. Indyk, and V. S. Mirrokni.  Locality-sensitive hashing scheme based on p-stable distributions.  In J. Snoeyink and J. Boissonnat, editors, Proceedings of the 20th ACM Symposium on Computational Geometry, Brooklyn, New York, USA, June 8-11, 2004, pages 253–262. ACM, 2004.
  * [10] A. Efimov, M. Gaal, and S. G. Revesz.  On integral estimates of non-negative positive definite functions, 2016.
  * [11] P. Flajolet and G. N. Martin.  Probabilistic counting algorithms for data base applications.  J. Comput. Syst. Sci., 31(2):182–209, sep 1985.
  * [12] M. G. Genton.  Classes of kernels for machine learning: A statistics perspective.  J. Mach. Learn. Res., 2:299–312, 2001.
  * [13] B. Giraud and R. Peschanski.  On positive functions with positive fourier transforms.  Acta Physica Polonica B, 37, 05 2005.
  * [14] B. Giraud and R. Peschanski.  From “dirac combs” to fourier-positivity.  Acta Physica Polonica B, 47(4):1075, 2016.
  * [15] B. G. Giraud and R. Peschanski.  On the positivity of fourier transforms, 2014.
  * [16] E. Hernández-Lemus.  Random fields in physics, biology and data science.  Frontiers in Physics, 9:77, 2021.
  * [17] A. Hinrichs and J. Vybíral.  On positive positive-definite functions and bochner’s theorem.  Journal of Complexity, 27(3):264–272, 2011.  Dagstuhl 2009.
  * [18] P. Indyk and R. Motwani.  Approximate nearest neighbors: Towards removing the curse of dimensionality.  In J. S. Vitter, editor, Proceedings of the Thirtieth Annual ACM Symposium on the Theory of Computing, Dallas, Texas, USA, May 23-26, 1998, pages 604–613. ACM, 1998.
  * [19] Z. Lin and Z. Bai.  Inequalities Related to Characteristic Functions, pages 23–28.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2011.
  * [20] C. Mercadier.  Numerical bounds for the distributions of the maxima of some one- and two-parameter gaussian processes.  Advances in Applied Probability, 38(1):149–170, 2006.
  * [21] O. R. Musin.  Multivariate positive definite functions on spheres.  Unkown, 1942.
  * [22] K. E. Pilario, M. Shafiee, Y. Cao, L. Lao, and S.-H. Yang.  A review of kernel methods for feature extraction in nonlinear process monitoring.  Processes, 8(1), 2020.
  * [23] C. E. Rasmussen and C. K. I. Williams.  Gaussian processes for machine learning.  Adaptive computation and machine learning. MIT Press, 2006.
  * [24] D. Slepian.  The one-sided barrier problem for gaussian noise.  The Bell System Technical Journal, 41(2):463–501, 1962.
  * [25] D. Slepian.  On the zeros of gaussian noise.  Time series analysis, pages 104–115, 1963.
  * [26] J. Verger-Gaugry.  Covering a ball with smaller equal balls in r${}^{\mbox{n}}$.  Discret. Comput. Geom., 33(1):143–155, 2005.
  * [27] G. N. Watson.  A treatise on the theory of bessel functions.  In A Treatise on the Theory of Bessel Functions. Cambridge University Press, Cambridge, 1966.