9
$\begingroup$

Let us sample some big number N points randomly uniformly on $[0,1]^d$. Consider 1-nearest neighbor graph based on such data cloud. (Let us look on it as UNdirected graph).

Question What would the number of connected components depending on d ? (As a percent of "N" - number of points.)

Simulation below suggest 31% for d=2, 20% for d=20, etc:

Percent Dimension:
31 2
28 5
25 10
20 20
15 50
13 100
10 1000

See code below. (One can run it on colab.research.google.com without installing anything on your comp).

If someone can comment on more general questions here: https://mathoverflow.net/q/362721/10446 that would be greatly appreciated.

!pip install python-igraph
!pip install cairocffi
import igraph

import time
from sklearn.neighbors import NearestNeighbors
import numpy as np

t0 = time.time()

dim = 20
n_sample = 10**4

for i in range(10): # repeat simulation 10 times to get more stat
  X = np.random.rand(n_sample, dim)
  nbrs = NearestNeighbors(n_neighbors=2, algorithm='ball_tree', ).fit(X)
  distances, indices = nbrs.kneighbors(X)
  g = igraph.Graph( directed = True )
  g.add_vertices(range(n_sample))
  g.add_edges(indices )
  g2 = g.as_undirected(mode = 'collapse')
  r = g2.clusters()
  print(len(r),len(r)/n_sample*100 , time.time() - t0)
$\endgroup$
  • 1
    $\begingroup$ FYI, for dimensions 10 and higher, the numbers reported above for the larger dimensions are inconsistent with the theoretical values calculated in my answer. Not sure where the discrepancy is coming from. Probably I've made a mistake, or you might need to increase n_sample for larger dimensions to approach the asymptotic result. $\endgroup$ – Neal Young Jun 14 at 15:09
16
$\begingroup$

For $n$ uniformly random points in a unit square the number of components is

$$\frac{3\pi}{8\pi+3\sqrt{3}}n+o(n)$$

See Theorem 2 of D. Eppstein, M. S. Paterson, and F. F. Yao (1997), "On nearest-neighbor graphs", Disc. Comput. Geom. 17: 263–282, https://www.ics.uci.edu/~eppstein/pubs/EppPatYao-DCG-97.pdf

For points in any fixed higher dimension it is $\Theta(n)$; I don't know the exact constant of proportionality but the paper describes how to calculate it.

| cite | improve this answer | |
$\endgroup$
  • 2
    $\begingroup$ The Theta is an instance of big-theta notation, a close relative of big-O. It means that it is both lower-bounded and upper-bounded by constant multiples of n. (In this case it's the same constant, but the notation implies only that there exist both lower-bound and upper-bound constants.) $\endgroup$ – David Eppstein Jun 13 at 16:32
  • 3
    $\begingroup$ For $k$NN graphs there are still linearly many components; the more interesting question is whether there is an infinite component. This turns out to be true for sufficiently large $k$ depending on the dimension; see Shang-Hua Teng and Frances F. Yao (2007), "$k$-Nearest-Neighbor Clustering and Percolation Theory", Algorithmica 49: 192–211, doi.org/10.1007/s00453-007-9040-7 $\endgroup$ – David Eppstein Jun 14 at 5:20
  • 1
    $\begingroup$ I don't know a reference for that version but I expect the answer to be $\Theta(\log n)$. There's a reasonable probability of getting a region of logarithmic area, containing a logarithmic number of points, surrounded by an empty region wider than the diameter of the inner region, so it's definitely at least $\Omega(\log n)$. $\endgroup$ – David Eppstein Jun 15 at 6:59
  • 1
    $\begingroup$ There aren't many papers citing the Teng-Yao paper so it's not hard to look through them and find what's known. There's a claimed improvement of k_2 to 188 (not 3) in Bagchi and Bansal, "Nearest-neighbor graphs on random point sets and their applications to sensor networks", PODC 2008, but without detail and with no subsequent detailed publication. I didn't see anything else relevant. $\endgroup$ – David Eppstein Jun 18 at 19:57
  • 1
    $\begingroup$ The second-order term is nonexistent for Poisson models in the whole plane, and stems only from the boundary effects in cutting off that plane to make a square and from the small differences between Poisson and sampling an exact number of points. $\endgroup$ – David Eppstein Jun 21 at 23:55
10
$\begingroup$

EDIT 2: Made explicit the underlying non-asymptotic bounds in the calculation.

EDIT: Replaced the calculation for two dimensions by the case of arbitrary constant dimension. Added a table of the values.

I'd like to add an informal sketch of how David's very elegant result can be calculated. (To be clear, I suggest selecting his answer as the "correct" answer; this one is just intended to complement his.)

Assume the points are in general position, so that no two distinct pairs have the same distance. This happens with probability 1.

  1. In the directed nearest-neighbor graph, each point has out-degree 1 (by definition). Also, for any directed path $p_1 \rightarrow p_2 \rightarrow p_3 \rightarrow \cdots \rightarrow p_k$ with no 2-cycles, we have $d(p_1, p_2) > d(p_2, p_3) > \cdots > d(p_{k-1}, p_k)$. That is, edge lengths decrease along the path. This is because, e.g., $p_3$ must be closer to $p_2$ than $p_1$ is (otherwise $p_3$ would not be $p_2$'s nearest neighbor), and so on.

  2. As a consequence, in the undirected multigraph obtained by replacing each directed edge by its undirected equivalent, the only cycles are 2-cycles, where points $p_i$ and $p_j$ form a 2-cycle if and only if they are mutual nearest neighbors. Other edges are not on cycles.

  3. It follows that the undirected nearest-neighbor graph (in which each such 2-cycle is replaced by a single edge) is acyclic, and has a number of edges equal to the number of vertices minus the number of pairs that are mutual nearest neighbors. Thus the number of components is equal to the number of pairs that are mutual nearest neighbors.

This holds in any metric space. Next, for intuition, consider the case of points in $R^1$, where calculations are relatively easy.

One Dimension

To ease the calculation, modify the distance metric to "wrap around" the boundary, that is, use $$d_1(x, x') = \min\{|x-x'|,1-|x-x'|\}.$$ This changes the number of nearest-neighbor pairs by at most 1.

We need to estimate the expected number of pairs, among the $n$ points, that are mutual nearest neighbors. If we order the points as $$p_1 < p_2 < \cdots < p_{n},$$ only pairs of the form $(p_i, p_{i+1})$ (or $(p_n, p_1)$) can be nearest neighbors. A given pair of this form are nearest-neighbors if and only if their distance is less than the distances of the neighboring pairs $(p_{i-1}, p_i)$ and $(p_{i+1}, p_{i+2})$ (to the left and right). That is, if, among the three consecutive pairs, the middle pair is closest. By symmetry (?), this happens with probability 1/3. Hence, by linearity of expectation, the number of the $n$ adjacent pairs that are nearest-neighbor pairs is $n/3$ (plus or minus 1, to correct for the wrap-around assumption). Hence the number of components is $n/3\pm 1$.

The symmetry argument above is suspect -- maybe there is some conditioning? It also doesn't extend to higher dimensions. Here's a more careful, detailed calculation that addresses these issues. Let $p_1, p_2, \ldots, p_n$ be the points in sampled order. By linearity of expectation, the expected number of nearest-neighbor pairs is the number of pairs $n\choose 2$ times the probability that a given random pair $(p, q)$ are a nearest-neighbor pair. WLOG we can assume that $p$ and $q$ are the first two points drawn. Let $d_{pq}$ be their distance. They will be a nearest-neighbor pair if and only if none of the $n-2$ subsequent points are within distance $d_{pq}$ of $p$ or $q$.

The probability of this event (conditioned on $d_{pq}$) is $$\max(0, 1-3d_{pq})^{n-2},$$ because it happens if and only if none of the remaining $n-2$ points falls between $p$ and $q$ or within either of the two $d_{pq}$-wide boundaries on each side of $p$ and $q$.

To distance $d_{pq}$ is distributed uniformly over $[0, 1/2]$ (using our "wrap-around" assumption). Hence, the probability that $(p,q)$ is a nearest-neighbor pair is $$\int_{0}^{1/3} (1-3 z)^{n-2} 2dz.$$ By a change of variables $x = 1-3z$ this is $$\int_{0}^1 x^{n-2} 2\,dx/3 = \frac{2}{3(n-1)}.$$ By linearity of expectation, the expected number of nearest-neighbor pairs is $2{n\choose 2}/(3(n-1)) = n/3$ (plus or minus 1 to correct for the wrap-around technicality). Therefore the expected number of components is indeed $n/3\pm 1$.

As an aside, note that when $d_{pq}$ is large (larger than $\log(n)/n$, say), the contribution to the expectation above is negligible. So, we could under- or over-estimate the conditional probability for such $d_{pq}$ significantly; this would change the result by lower-order terms.

Any constant dimension

Fix constant dimension $k \in \{1,2,\ldots\}$.

To ease the calculations, modify the distance metric to wrap around the borders, that is, use $d_k(p, q) = \sqrt{\sum_{i=1}^k d_1(p_i, q_i)^2}$ for $d_1$ as defined previously. This changes the answer by at most an additive $o(n)$ with high probability and in expectation.

Define $\beta_k, \mu_k\in \mathbb R$ such that $\beta_k r^k$ and $\mu_k r^k$ are the volumes of, respectively, a ball of radius $r$ and the union of two overlapping balls of radius $r$ whose centers are $r$ apart (so each center lies on the boundary of the other ball).

Let $p_1, p_2, \ldots, p_n$ be the points in sampled order. By linearity of expectation, the expected number of nearest-neighbor pairs is the number of pairs $n\choose 2$ times the probability that a given random pair $(p, q)$ are a nearest-neighbor pair. WLOG we can assume that $p$ and $q$ are the first two points drawn. Let $d_{pq}$ be their distance. They will be a nearest-neighbor pair if and only if none of the $n-2$ subsequent points are within distance $d_{pq}$ of $p$ or $q$.

We calculate the probability of this event. In the case that $d_{pq} \ge 1/4$, the probability of the event is at most the probability that no point falls within the ball of radius $1/4$ around $p$, which is at most $(1-\beta_k/4^k)^{n-2} \le \exp(-(n-2)\beta_k/4^k)$.

The case that $d_{pq} \le 1/4$ happens with probability $\beta_k/4^k$. Condition on any such $d_{pq}$. Then $p$ and $q$ will be nearest neighbors iff none of the $n-2$ subsequent points fall in the "forbidden" region consisting of the union of the two balls of radius $d_{pq}$ with centers at $p$ and $q$. The area of this region is $\mu_k d_{pq}^k$ by definition of $\mu_k$ (using here that $d_{pq}\le 1/4$ and the metric wraps around), so the probability of the event in question is $(1-\mu_k d_{pq}^k)^{n-2}$.

Conditioned on $d_{pq} \in [0,1/4]$, the probability density function of $d_{pq}$ is $f(r) = k 4^k r^{k-1}$ (note $\int_{0}^{1/4} k 4^k r^{k-1} = 1$). Hence, the overall (unconditioned) probability of the event is $$\frac{\beta_k}{4^k} \int_{0}^{1/4} (1-\mu_k r^k)^{n-2} k 4^kr^{k-1} \, dr ~+~ \epsilon(n,k)$$ where $$0 \le \epsilon(n, k) \le \exp(-(n-2)\beta_k /4^k).$$ Using a change of variables $z^k=1-\mu_k r^k$ to calculate the integral, this is $$\frac{k \beta_k}{\mu_k} \int_{\alpha}^1 z^{k(n-1)-1} \, dz ~+~ \epsilon'(n, k) = \frac{\beta_k}{\mu_k}\frac{1 + \epsilon'(n, k)}{n-1}$$ for constant $\alpha=(1-\mu_k/4^k)^{1/k}<1$ and "error term" $\epsilon'(n, k)$ satisfying $$-\exp(-(n-1)\mu_k/4^k) ~\le~ \epsilon'(n, k) ~\le~ \exp(-(n-2)\beta_k /4^k)(n-1)\mu_k/\beta_k$$ so $\epsilon'(n, k) \rightarrow 0$ as $n\rightarrow\infty$.

By linearity of expectation, the expected number of nearest-neighbor pairs is $$\frac{n\choose 2}{n-1}\frac{\beta_k}{\mu_k}(1+ \epsilon'(n,k)) = \frac{\beta_k}{2\mu_k}(1 + \epsilon'(n,k)) n,$$ where $\epsilon'(n, k) \rightarrow 0$ as $n\rightarrow\infty$. Correcting for the wrap-around assumption adds a $\pm o(n)$ term.

Hence, asymptotically, the expected number of mutual nearest-neighbor pairs is $n\beta_k/(2\mu_k) + o(n)$. Next we give more explicit forms for $\beta_k/(2\mu_k)$.

According to this Wikipedia entry, $$\beta_k = \frac{\pi^{k/2}}{\Gamma(k/2+1)} \sim \frac{1}{\sqrt{\pi k}}\Big(\frac{2\pi e}{k}\Big)^{k/2}$$ where $\Gamma$ is Euler's gamma function, with $\Gamma(k/2+1) \sim \sqrt{\pi k}(k/(2e))^{k/2}$ (see here).

Following the definition of $\mu_k$, the volume of the union of the two balls is twice the volume of one ball with a "cap" removed (where the cap contains the points in the ball that are closer to the other ball). Using this math.se answer (taking $d=r_1=r_2=r$, so $c_1=a=r/2$) to get the volume of the cap, this gives $$\mu_k = \beta_k (2 - I_{3/4}((k+1)/2, 1/2)),$$ where $I$ is the "regularized incomplete beta function". Hence, the desired ratio is $$\frac{\beta_k}{2\mu_k} = \Big(4-2I_{3/4}\Big(\frac{k+1}2, \frac{1}{2}\Big)\Big)^{-1}.$$

I've appended below the first 20 values, according to WolframAlpha.


$$ \begin{array}{cc} \begin{array}{|rcl|} k & \beta_k / (2\mu_k) & \approx \\ \hline 1 & \displaystyle\frac{1}{3} & 0.333333 \\ 2 &\displaystyle\frac{3 \pi}{3 \sqrt{3}+8 \pi} & 0.310752 \\ 3 &\displaystyle\frac{8}{27} & 0.296296 \\ 4 &\displaystyle\frac{6 \pi}{9 \sqrt{3}+16 \pi} & 0.286233 \\ 5 &\displaystyle\frac{128}{459} & 0.278867 \\ 6 &\displaystyle\frac{15 \pi}{27 \sqrt{3}+40 \pi} & 0.273294 \\ 7 &\displaystyle\frac{1024}{3807} & 0.268978 \\ 8 &\displaystyle\frac{420 \pi}{837 \sqrt{3}+1120 \pi} & 0.265577 \\ 9 &\displaystyle\frac{32768}{124659} & 0.262861 \\ 10 &\displaystyle\frac{420 \pi}{891 \sqrt{3}+1120 \pi} & 0.26067 \\ \hline \end{array} & \begin{array}{|rcl|} k & \beta_k / (2\mu_k) & \approx \\ \hline 11 &\displaystyle\frac{262144}{1012581} & 0.258887 \\ 12 &\displaystyle\frac{330 \pi}{729 \sqrt{3}+880 \pi} & 0.257427 \\ 13 &\displaystyle\frac{4194304}{16369695} & 0.256224 \\ 14 &\displaystyle\frac{5460 \pi}{12393 \sqrt{3}+14560 \pi} & 0.255228 \\ 15 &\displaystyle\frac{33554432}{131895783} & 0.254401 \\ 16 &\displaystyle\frac{120120 \pi}{277749 \sqrt{3}+320320 \pi} & 0.253712 \\ 17 &\displaystyle\frac{2147483648}{8483550147} & 0.253135 \\ 18 &\displaystyle\frac{2042040 \pi}{4782969 \sqrt{3}+5445440 \pi} & 0.252652 \\ 19 &\displaystyle\frac{17179869184}{68107648041} & 0.252246 \\ 20 & \displaystyle\frac{38798760 \pi}{91703097 \sqrt{3}+103463360 \pi} & 0.251904 \\ \hline \end{array} \end{array} $$


Notably, for $k=20$ (and larger, in fact), WolframAlpha reports numerical values approaching 0.25, in contrast to the experimental results reported by OP which are much lower. Where is this discrepancy coming from?

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ We see that percent is growing with point number. But still for D=20 the gap to theory is big. So it is quite interesting puzzle. $\endgroup$ – Alexander Chervov Jun 19 at 10:31
  • 1
    $\begingroup$ I wonder if the calculation in my answer has an error.. I checked again but still haven't found any error. Running your script on colab.google, with $n=10^5$, the observed value is slightly smaller than the theoretical value starting with dimension $k=8$ or so, it's about 1% low (and less than 25%) already for dimension $k=10$. But it's not clear that the asymptotically small error terms in the calculation are actually small when $k=10$ and $n=10^6$.. I think a next step would be to more carefully consider the error terms in the calculation, especially the negative ones. $\endgroup$ – Neal Young Jul 7 at 14:18
  • 1
    $\begingroup$ @AlexanderChervov, I edited my answer to give non-asymptotic upper and lower bounds on the (asymptotically vanishing) error terms. The lower bounds do not rule out a constant-size gap (error term), as you may be seeing in your calculations, until $n$ is over the threshold $$n \gg 4^k/\beta_k = (8/\pi)^{k/2} \Gamma(k/2+1).$$ This threshold is over 1M for $k=14$, according to Wolfram Alpha. But perhaps refining the analysis could give tighter bounds. $\endgroup$ – Neal Young Jul 8 at 2:45
  • 1
    $\begingroup$ @AlexanderChervov, I tried some (possibly naive) way of bounding the error introduced by assuming the metric wraps around, which I hadn't considered in my previous comment. The bound I got does not rule out a constant-size gap until $n$ is over the threshold $n \gg k^k/\beta_k$, even larger than the threshold mentioned in my previous comment. (Which, FWIW, I think I can reduce from $4^k/\beta_k$ to $2^k/\beta_k$. So the remaining bottleneck is bounding the error introduced by the "wrap-around" assumption.) $\endgroup$ – Neal Young Jul 8 at 15:27
  • 1
    $\begingroup$ It might be possible to estimate the "wrap-around" error empirically by redoing some of your experiments but changing the distance metric to $d_k$ as defined in my answer. If you do that, I think you should observe that the gap is small for $n\gg 2^k/\beta_k \sim \sqrt{\pi k}(\frac{k}{4.27})^{k/2}$. $\endgroup$ – Neal Young Jul 8 at 15:34
0
$\begingroup$

It is not an answer, but a comment to very interesting answer by Neal Young; who proposed the beautiful formula for general dimensions; thus generalizing the beautiful formula by David Eppstein and coauthors. The formula fits simulations very nicely for lower dimensions; however the descrepancy appears for higher dimensions.

Thus it is quite interesting to understand the reason for the descrepancy. It might have interesting applied consequences, for example for testing KNN algorithms and their approximate versions.

There, of course, can be simple reasons, for the discrepancy - some of us made a mistake, but can be also more interesting reason - like - I am not simulating large enough number of points (but see below even for 10 millions) - thus interesting to understand next term in asymptotic, or some other.

So let me share some more simulations results, self-checks, comments and scripts.

In conclusion: it seems simulations are correct, at least I checked several issues (actually not all of them) I have worried about. For large dimensions like 50 we have quite a big discrepancy, it would be very interesting if it can be explained by a small sample size, that would imply existence of extremely powerful second order term...

Actually simulating big sizes is somewhat tricky, I am still not 100 percent sure in correctness. Probably the main point of writing all these - is to share possible subtleties which might appear, if someone repeats the simulations.

Dimension = 10 , Theory Percent 26.067

Sample Size   Percent by Simulation
1 000         24.1311
10 000        24.5819
100 000       24.90933
1 000 000     25.146969
10 000 000    25.342639

We see result slightly grows with sample size, (however for large dimensions it would not be true) So it might be that increasing size we get agreement with theory, though growth is quite small. The simulation is done repeatedly 100 times (except last size where only 10 times). The script can be found here: https://www.kaggle.com/alexandervc/connected-components-knn-graph-v010-rapids-knn?scriptVersionId=38115858 The simulation is using GPU package RAPIDs based on Facebooks FAISS https://engineering.fb.com/data-infrastructure/faiss-a-library-for-efficient-similarity-search/ GPU can accelerate these calculations up to 500 times. It is done on kaggle platform, where you can use 9 nine hours of GPU continuously and 30 hours as whole per week for free, and where all these GPU packages are installable correctly. Many thanks to Dmitry Simakov for sharing his notebooks, letting me know about RAPIDs etc.

What is subtle here: It is known that GPU is single precision, while CPU is double precesion - and surprisingly enough that causes small difference in graphs produced. (It is known). However this small numerical instability should not affect statistical properties. I hope so , or it might be interesting point that it is not like that.

Dimension = 50 , Theory Percent about 25

Sample Size   Percent by Simulation
1 000         16.694
10 000        15.6265
100 000       15.05882
1 000 000     14.834492

Notebook: https://www.kaggle.com/alexandervc/connected-components-knn-graph-v010-rapids-knn?scriptVersionId=38115858

We see that even increasing the sample size the percent does not increase and it is quite far from theory. Again see subtlety mentioned above.

What is subtle here: See above

Dimension = 20 , Theory Percent about 25.19

Sample Size   Percent by Simulation
1 000         21.3
10 000        20.15
100 000       20.817
1 000 000     21.3472
10 000 000    21.817

There is small increase with sample size, but theory is quite far...

Notebook up to 1 000 000 : https://www.kaggle.com/alexandervc/connected-components-knn-graph-v010-rapids-knn?scriptVersionId=37225738 Notebook for 10 000 000 : https://www.kaggle.com/alexandervc/connected-components-knn-graph-v010-rapids-knn?scriptVersionId=37148875

Dimension = 5 (100 times averaged) Percent Theory = 27.8867

Size Mean         Std/sqrt(100)
1e3 27.531000 +- 0.0720787
1e4 27.650000 +- 0.0255797
1e5 27.745390 +- 0.0069290 
1e6 27.794086 +- 0.0024427
1e7 27.830324 +- 0.00072 

1e7 - time: 446.144 seconds - per 1 run 1e6 - time: 26.1098 seconds - per 1 run

What is subtle here: That simulation is done on colab CPU, the point is one can use NOT the brute force method to calculate KNN, graph, but kd_tree method (built-in in Python sklearn), which is exact (not approximate), but works much faster than brute force method which is scales quadratically with sample size. The problem is that it works fast for low dimensions like 5 (for uniform data), and begin to work MUCH SLOWER for higher dimensions.

Here is notebook with speed comparaison: https://www.kaggle.com/alexandervc/compare-nn-graph-speed-sklearn-vs-gpu-rapids

PS

I also checked the calculation of connected component count implemented by different Python packages - igraph and snap and networkX gives the same result. So it should not be a error at that part.

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ @NealYoung some new observation: consider sphere or torus - then results seems to agree with your proposal ! For example - sphere - dimension 10 , simulations up 10 millions - get 26.32 - a little ABOVE theory (not below as was before) - kaggle.com/alexandervc/… - please look tables at the end of notebook , I will report more later ... $\endgroup$ – Alexander Chervov Jul 11 at 20:19
  • 1
    $\begingroup$ @NealYoung Sphere Dimension 50 , gives coefficient 25.05 - which finely agrees, with you prediction that results should be above 25 for any dimension. kaggle.com/alexandervc/… $\endgroup$ – Alexander Chervov Jul 12 at 7:55
  • 1
    $\begingroup$ That's a relief. So it seems that for the hypercube it is the terms due to boundary effects that vanish most slowly as $n$ grows. It would be an interesting exercise to try to quantify that more precisely. $\endgroup$ – Neal Young Jul 13 at 2:29
  • 1
    $\begingroup$ @NealYoung sphere dimensions 500 and 5000 also give coefficent around 25 within std accuracy. For 500 kaggle.com/alexandervc/… and for 5000 kaggle.com/alexandervc/… $\endgroup$ – Alexander Chervov Jul 13 at 8:09
  • $\begingroup$ Even dimension 50 000 is also within simulation possibility (for sample sizes like 1000-5000 , and may be more) and again results are around 25 percent in full agreement with the theory. kaggle.com/alexandervc/… $\endgroup$ – Alexander Chervov Jul 14 at 17:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.