<p><em>EDIT 2: Made explicit the underlying non-asymptotic bounds in the calculation.</em></p>
<p><em>EDIT: Replaced the calculation for two dimensions by the case of arbitrary constant dimension.  Added a table of the values.</em></p>
<p>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 &quot;correct&quot; answer; this one is just intended to complement his.)</p>
<p>Assume the points are in general position, so that no two distinct pairs have the same distance.  This happens with probability 1.</p>
<ol>
<li><p>In the directed nearest-neighbor graph, each point has out-degree 1 (by definition).  Also, for any directed path <span class="math-container">$p_1 \rightarrow p_2 \rightarrow p_3 \rightarrow \cdots \rightarrow  p_k$</span> with no 2-cycles, we have <span class="math-container">$d(p_1, p_2) &gt; d(p_2, p_3) &gt; \cdots &gt; d(p_{k-1}, p_k)$</span>.  That is, edge lengths decrease along the path.  This is because, e.g., <span class="math-container">$p_3$</span> must be closer to <span class="math-container">$p_2$</span> than <span class="math-container">$p_1$</span> is (otherwise <span class="math-container">$p_3$</span> would not be <span class="math-container">$p_2$</span>'s nearest neighbor), and so on.</p>
</li>
<li><p>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 <span class="math-container">$p_i$</span> and <span class="math-container">$p_j$</span> form a 2-cycle if and only if they are mutual nearest neighbors.  Other edges are not on cycles.</p>
</li>
<li><p>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 <em>the number of components is equal to the number of pairs that are mutual nearest neighbors.</em></p>
</li>
</ol>
<p>This holds in any metric space.  Next, for intuition, consider the case of points in <span class="math-container">$R^1$</span>, where calculations are relatively easy.</p>
<h2>One Dimension</h2>
<p>To ease the calculation, modify the distance metric to &quot;wrap around&quot; the boundary, that is, use
<span class="math-container">$$d_1(x, x') = \min\{|x-x'|,1-|x-x'|\}.$$</span>  This changes the number of nearest-neighbor pairs by at most 1.</p>
<p>We need to estimate the expected number of pairs, among the <span class="math-container">$n$</span> points, that are mutual nearest neighbors. If we order the points as
<span class="math-container">$$p_1 &lt; p_2 &lt; \cdots &lt; p_{n},$$</span>
only pairs of the form <span class="math-container">$(p_i, p_{i+1})$</span> (or <span class="math-container">$(p_n, p_1)$</span>) 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 <span class="math-container">$(p_{i-1}, p_i)$</span> and <span class="math-container">$(p_{i+1}, p_{i+2})$</span> (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 <span class="math-container">$n$</span> adjacent pairs that are nearest-neighbor pairs is <span class="math-container">$n/3$</span> (plus or minus 1, to correct for the wrap-around assumption). Hence the number of components is <span class="math-container">$n/3\pm 1$</span>.</p>
<p>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 <span class="math-container">$p_1, p_2, \ldots, p_n$</span> be the points in sampled order. By linearity of expectation, the expected number of nearest-neighbor pairs is the number of pairs <span class="math-container">$n\choose 2$</span> times the probability that a given random pair <span class="math-container">$(p, q)$</span> are a nearest-neighbor pair.  WLOG we can assume that <span class="math-container">$p$</span> and <span class="math-container">$q$</span> are the first two points drawn.  Let <span class="math-container">$d_{pq}$</span> be their distance.  They will be a nearest-neighbor pair if and only if none of the <span class="math-container">$n-2$</span> subsequent points are within distance <span class="math-container">$d_{pq}$</span> of <span class="math-container">$p$</span> or <span class="math-container">$q$</span>.</p>
<p>The probability of this event (conditioned on <span class="math-container">$d_{pq}$</span>) is <span class="math-container">$$\max(0, 1-3d_{pq})^{n-2},$$</span>
because it happens if and only if none of the remaining <span class="math-container">$n-2$</span> points falls between <span class="math-container">$p$</span> and <span class="math-container">$q$</span> or within either of the two <span class="math-container">$d_{pq}$</span>-wide boundaries on each side of <span class="math-container">$p$</span> and <span class="math-container">$q$</span>.</p>
<p>To distance <span class="math-container">$d_{pq}$</span> is distributed uniformly over <span class="math-container">$[0, 1/2]$</span> (using our &quot;wrap-around&quot; assumption).  Hence, the probability that <span class="math-container">$(p,q)$</span> is a nearest-neighbor pair is
<span class="math-container">$$\int_{0}^{1/3} (1-3 z)^{n-2} 2dz.$$</span>
By a change of variables <span class="math-container">$x = 1-3z$</span> this is
<span class="math-container">$$\int_{0}^1 x^{n-2} 2\,dx/3 = \frac{2}{3(n-1)}.$$</span>
By linearity of expectation, the expected number of nearest-neighbor pairs is <span class="math-container">$2{n\choose 2}/(3(n-1)) = n/3$</span> (plus or minus 1 to correct for the wrap-around technicality).  Therefore the expected number of components is indeed <span class="math-container">$n/3\pm 1$</span>.</p>
<p>As an aside, note that when <span class="math-container">$d_{pq}$</span> is large (larger than <span class="math-container">$\log(n)/n$</span>, say), the contribution to the expectation above is negligible.  So, we could under- or over-estimate the conditional probability for such <span class="math-container">$d_{pq}$</span> significantly; this would change the result by lower-order terms.</p>
<h2>Any constant dimension</h2>
<p>Fix constant dimension <span class="math-container">$k \in \{1,2,\ldots\}$</span>.</p>
<p>To ease the calculations, modify the distance metric to wrap around the borders, that is, use
<span class="math-container">$d_k(p, q) = \sqrt{\sum_{i=1}^k d_1(p_i, q_i)^2}$</span>
for <span class="math-container">$d_1$</span> as defined previously.  This changes the answer by at most an additive <span class="math-container">$o(n)$</span> with high probability and in expectation.</p>
<p>Define <span class="math-container">$\beta_k, \mu_k\in \mathbb R$</span> such that <span class="math-container">$\beta_k r^k$</span> and <span class="math-container">$\mu_k r^k$</span> are the volumes of, respectively, a ball of radius <span class="math-container">$r$</span> and the union of two overlapping balls of radius <span class="math-container">$r$</span> whose centers are <span class="math-container">$r$</span> apart (so each center lies on the boundary of the other ball).</p>
<p>Let <span class="math-container">$p_1, p_2, \ldots, p_n$</span> be the points in sampled order. By linearity of expectation, the expected number of nearest-neighbor pairs is the number of pairs <span class="math-container">$n\choose 2$</span> times the probability that a given random pair <span class="math-container">$(p, q)$</span> are a nearest-neighbor pair.  WLOG we can assume that <span class="math-container">$p$</span> and <span class="math-container">$q$</span> are the first two points drawn.  Let <span class="math-container">$d_{pq}$</span> be their distance.  They will be a nearest-neighbor pair if and only if none of the <span class="math-container">$n-2$</span> subsequent points are within distance <span class="math-container">$d_{pq}$</span> of <span class="math-container">$p$</span> or <span class="math-container">$q$</span>.</p>
<p>We calculate the probability of this event. In the case that <span class="math-container">$d_{pq} \ge 1/4$</span>, the probability of the event is at most the probability that no point falls within the ball of radius <span class="math-container">$1/4$</span> around <span class="math-container">$p$</span>, which is at most <span class="math-container">$(1-\beta_k/4^k)^{n-2} \le \exp(-(n-2)\beta_k/4^k)$</span>.</p>
<p>The case that <span class="math-container">$d_{pq} \le 1/4$</span> happens with probability <span class="math-container">$\beta_k/4^k$</span>.  Condition on any such <span class="math-container">$d_{pq}$</span>.  Then <span class="math-container">$p$</span> and <span class="math-container">$q$</span> will be nearest neighbors iff none of the <span class="math-container">$n-2$</span> subsequent points fall in the &quot;forbidden&quot; region consisting of the union of the two balls of radius <span class="math-container">$d_{pq}$</span> with centers at <span class="math-container">$p$</span> and <span class="math-container">$q$</span>.  The area of this region is <span class="math-container">$\mu_k d_{pq}^k$</span> by definition of <span class="math-container">$\mu_k$</span> (using here that <span class="math-container">$d_{pq}\le 1/4$</span> and the metric wraps around), so the probability of the event in question is <span class="math-container">$(1-\mu_k d_{pq}^k)^{n-2}$</span>.</p>
<p>Conditioned on <span class="math-container">$d_{pq} \in [0,1/4]$</span>, the probability density function of <span class="math-container">$d_{pq}$</span> is <span class="math-container">$f(r) = k 4^k r^{k-1}$</span> (note <span class="math-container">$\int_{0}^{1/4} k 4^k r^{k-1} = 1$</span>). Hence, the overall (unconditioned) probability of the event is
<span class="math-container">$$\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)$$</span>
where <span class="math-container">$$0 \le \epsilon(n, k) \le \exp(-(n-2)\beta_k /4^k).$$</span>
Using a change of variables <span class="math-container">$z^k=1-\mu_k r^k$</span> to calculate the integral, this is
<span class="math-container">$$\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}$$</span>
for constant <span class="math-container">$\alpha=(1-\mu_k/4^k)^{1/k}&lt;1$</span>
and &quot;error term&quot; <span class="math-container">$\epsilon'(n, k)$</span> satisfying
<span class="math-container">$$-\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$$</span>
so <span class="math-container">$\epsilon'(n, k) \rightarrow 0$</span> as <span class="math-container">$n\rightarrow\infty$</span>.</p>
<p>By linearity of expectation, the expected number of nearest-neighbor pairs is
<span class="math-container">$$\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,$$</span>
where <span class="math-container">$\epsilon'(n, k) \rightarrow 0$</span> as <span class="math-container">$n\rightarrow\infty$</span>.
Correcting for the wrap-around assumption adds a <span class="math-container">$\pm o(n)$</span> term.</p>
<p>Hence, asymptotically, the expected number of mutual nearest-neighbor pairs is <span class="math-container">$n\beta_k/(2\mu_k) + o(n)$</span>.
Next we give more explicit forms for <span class="math-container">$\beta_k/(2\mu_k)$</span>.</p>
<p>According to <a href="https://en.wikipedia.org/wiki/Volume_of_an_n-ball#The_volume" rel="nofollow noreferrer">this</a> Wikipedia entry,
<span class="math-container">$$\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}$$</span>
where <span class="math-container">$\Gamma$</span> is Euler's gamma function,
with <span class="math-container">$\Gamma(k/2+1) \sim \sqrt{\pi k}(k/(2e))^{k/2}$</span> (see <a href="https://en.wikipedia.org/wiki/Gamma_function#Stirling%27s_formula" rel="nofollow noreferrer">here</a>).</p>
<p>Following the definition of <span class="math-container">$\mu_k$</span>, the volume of the union of the two balls is twice the volume of one ball with a &quot;cap&quot; removed (where the cap contains the points in the ball that are closer to the other ball).
Using <a href="https://math.stackexchange.com/a/162873/48535">this</a> math.se answer (taking <span class="math-container">$d=r_1=r_2=r$</span>, so <span class="math-container">$c_1=a=r/2$</span>) to get the volume of the cap, this gives
<span class="math-container">$$\mu_k = \beta_k (2 - I_{3/4}((k+1)/2, 1/2)),$$</span>
where <span class="math-container">$I$</span> is the &quot;regularized incomplete beta function&quot;.  Hence, the desired ratio is
<span class="math-container">$$\frac{\beta_k}{2\mu_k} = \Big(4-2I_{3/4}\Big(\frac{k+1}2, \frac{1}{2}\Big)\Big)^{-1}.$$</span></p>
<p>I've appended below the first 20 values, <a href="https://www.wolframalpha.com/input/?i=1+%2F%284-+2+BetaRegularized%5B3%2F4%2C+%28k%2B1%29%2F2%2C+1%2F2%5D%29+for+k%3D1%2C2%2C3%2C4%2C5%2C6%2C7%2C8%2C9%2C10%2C11%2C12%2C13%2C14%2C15%2C16%2C17%2C18%2C19%2C20" rel="nofollow noreferrer">according to WolframAlpha</a>.</p>
<hr />
<p><span class="math-container">$$
\begin{array}{cc}
\begin{array}{|rcl|}
k &amp; \beta_k / (2\mu_k) &amp; \approx \\ \hline
1 &amp; \displaystyle\frac{1}{3} &amp;  0.333333 \\
2 &amp;\displaystyle\frac{3 \pi}{3 \sqrt{3}+8 \pi} &amp;  0.310752 \\
3 &amp;\displaystyle\frac{8}{27} &amp;  0.296296 \\
4 &amp;\displaystyle\frac{6 \pi}{9 \sqrt{3}+16 \pi} &amp;  0.286233 \\
5 &amp;\displaystyle\frac{128}{459} &amp;  0.278867 \\
6 &amp;\displaystyle\frac{15 \pi}{27 \sqrt{3}+40 \pi} &amp;  0.273294 \\
7 &amp;\displaystyle\frac{1024}{3807} &amp;  0.268978 \\
8 &amp;\displaystyle\frac{420 \pi}{837 \sqrt{3}+1120 \pi} &amp;  0.265577 \\
9 &amp;\displaystyle\frac{32768}{124659} &amp;  0.262861 \\
10 &amp;\displaystyle\frac{420 \pi}{891 \sqrt{3}+1120 \pi} &amp;  0.26067 \\ \hline
\end{array} &amp;
\begin{array}{|rcl|}
k &amp; \beta_k / (2\mu_k) &amp; \approx \\ \hline
11 &amp;\displaystyle\frac{262144}{1012581} &amp;  0.258887 \\
12 &amp;\displaystyle\frac{330 \pi}{729 \sqrt{3}+880 \pi} &amp;  0.257427 \\
13 &amp;\displaystyle\frac{4194304}{16369695} &amp;  0.256224 \\
14 &amp;\displaystyle\frac{5460 \pi}{12393 \sqrt{3}+14560 \pi} &amp;  0.255228 \\
15 &amp;\displaystyle\frac{33554432}{131895783} &amp;  0.254401 \\
16 &amp;\displaystyle\frac{120120 \pi}{277749 \sqrt{3}+320320 \pi} &amp;  0.253712 \\
17 &amp;\displaystyle\frac{2147483648}{8483550147} &amp;  0.253135 \\
18 &amp;\displaystyle\frac{2042040 \pi}{4782969 \sqrt{3}+5445440 \pi} &amp;  0.252652 \\
19 &amp;\displaystyle\frac{17179869184}{68107648041} &amp;  0.252246 \\
20 &amp; \displaystyle\frac{38798760 \pi}{91703097 \sqrt{3}+103463360 \pi} &amp;  0.251904 \\ \hline
\end{array}
\end{array}
$$</span></p>
<hr />
<p>Notably, for <span class="math-container">$k=20$</span> (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?</p>

CC BY-SA 4.0
2020-06-14T04:14:08.267Z
very insightful answer, thank you :)
Sandeep Silwal
CC BY-SA 4.0
2020-06-14T16:04:22.867Z
Thank you very much ! I will double check my simulation
Alexander Chervov
CC BY-SA 4.0
2020-06-14T23:46:32.237Z
@AlexanderChervov: If the analysis above is right, the number has to be between 1/4 and 1/2 for every dimension because $\mu_k$ is the volume of the union of two intersecting balls of volume $\beta_k$ each (so $\beta_k/(2\mu_k)$ has to be between 1/4 and 1/2).  So either the analysis or your code is wrong.  Re: the code, are you sure it's using the $\ell_2$ metric?  Could precision issues somehow be a problem?  One sanity check would be to check that for every vertex $u$ in the directed graph, the path $u\rightarrow v \rightarrow w$ (of nearest neighbors) has $d(v, w) \le d(v, u)$.
Neal Young
CC BY-SA 4.0
2020-06-15T00:45:00.903Z
Also I'm fairly sure that $\lim_{k\rightarrow\infty} \beta_k/(2\mu_k) = 1/4$, because the  the volume of the "cap" becomes insignificant relative to the full ball as the dimension grows, so $\mu_k \sim 2\beta_k$.
Neal Young
CC BY-SA 4.0
2020-06-15T02:41:41.590Z
FYI I sanity-checked your code and didn't find any mistake or obvious inconsistency.
Neal Young
CC BY-SA 4.0
2020-06-15T08:10:21.080Z
Thank you for your remarks. I made preliminary checks, changed igraph to snap , got the same result . A friend of mine, using GPU studied greater size, like 1e6-7 points, he sees some increase of the answers for dim = 20 , but it is very low. Latter we will check a bit more.
Alexander Chervov
CC BY-SA 4.0
2020-06-19T10:27:53.190Z
# Dim = 20 

# we are quite far from theory - 25.1904
Alexander Chervov
CC BY-SA 4.0
2020-06-19T10:29:14.237Z
n: 10**4 points -  mean percent : 20.486, time: 19.92 s
Alexander Chervov
CC BY-SA 4.0
2020-06-19T10:29:35.353Z
n: 10**5, mean percent: 20.897, time: 60.168 s
Alexander Chervov
CC BY-SA 4.0
2020-06-19T10:30:06.777Z
n: 2*10**5 points , mean percent: 21.044, time: 148.198 s
Alexander Chervov
CC BY-SA 4.0
2020-06-19T10:31:15.267Z
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.
Alexander Chervov
CC BY-SA 4.0
2020-06-19T10:32:12.837Z
Later I will report on further simulations
Alexander Chervov
CC BY-SA 4.0
2020-07-07T06:45:16.717Z
Simulated dimension 50, it seems  percent is NOT growing with respect to sample size. n=100K 15.059 n=1Million 14.834 , with 100 times simulated and small std. Later I will report more. Need to check certain things
Alexander Chervov
CC BY-SA 4.0
2020-07-07T14:18:35.447Z
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.
Neal Young
CC BY-SA 4.0
2020-07-07T17:36:58.070Z
Thank you for your comment ! I have posted in answer my simulations with links and some comments (just in case). The main problem as it seems to me - for high dimensions - see D = 50 ( second example ). I will try to think more on that ....
Alexander Chervov
CC BY-SA 4.0
2020-07-08T02:45:52.297Z
@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](https://www.wolframalpha.com/input/?i=N%5B%288%2FPi%29%5E%28k%2F2%29+Gamma%28k%2F2%2B1%29%5D+where+k%3D14).  But perhaps refining the analysis could give tighter bounds.
Neal Young
CC BY-SA 4.0
2020-07-08T15:27:02.657Z
@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.)
Neal Young
CC BY-SA 4.0
2020-07-08T15:34:26.777Z
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}$.
Neal Young
CC BY-SA 4.0
2020-07-08T15:55:49.100Z
@NealYoung so how large should be "n" ? about ?
Alexander Chervov
CC BY-SA 4.0
2020-07-08T16:00:13.517Z
Python sklearn built-in functions should support several metrics: https://scikit-learn.org/stable/modules/generated/sklearn.neighbors.DistanceMetric.html?highlight=distancemetric%20class# But it seems you consider the torus (periodic) not the uniform cube , let me think how that can be done
Alexander Chervov
CC BY-SA 4.0
2020-07-08T16:32:48.397Z
Let us [continue this discussion in chat](https://chat.stackexchange.com/rooms/110359/discussion-between-neal-young-and-alexander-chervov).
Neal Young
CC BY-SA 4.0
2020-07-08T16:33:51.937Z
What I can do just use take points in (S^1 \in R^2)^dimension - that corresponds to yours metric, is it correct ?
Alexander Chervov
CC BY-SA 4.0
2020-12-03T09:23:20.150Z
If you have some time please take a look: https://cstheory.stackexchange.com/questions/47957/power-law-for-degree-distribution-of-random-knn-graphs
Alexander Chervov
CC BY-SA 4.0
2020-12-03T14:02:06.097Z
@AlexanderChervov I took a look and I don't know.. My guess would be that the maximum in-degree of any node is bounded by some function f(d, k) that depends only on d and k, independently of the number of points or how they are located (randomly, adversarially, ...).  Can you construct a counterexample to this guess?  Hopefully somebody who actually knows something about the literature (unlike me) will have something to say... :-)
Neal Young
CC BY-SA 4.0
2020-12-03T16:34:05.227Z
@NealYoung thank you very much for your answer I will try to check.
Alexander Chervov
CC BY-SA 4.0
2020-12-10T19:58:28.557Z
@NealYoung I tried to check whether maximum is bounded - it seems it is not (at least for Gaussian, dim =50) - see https://stats.stackexchange.com/questions/500250/recognize-distribution-on-integers-0-mean-5-median-2-percentiles-75-5   The distribution seems to have some nice properties like std apprximately 2*mean, mean = percentile 75
Alexander Chervov
CC BY-SA 4.0
2020-06-13T08:03:31.073Z
2020-07-08T13:16:13.827Z
2020-07-08T13:16:13.827Z
Neal Young
<p><strong>Update:</strong></p>
<p>Not a answer but sharing simulation result confirming proposed answer by Neal Young for certain cases.</p>
<p>Neal Young in his answer above a proposed beautiful formula for the number of connected components of NN-graph in any dimension. Simulation results in the table below quite confirm his proposal for the case of points distributed uniformly on sphere and torus (probably any closed manifold). Though for points taken from the uniform (on cube) or Gaussian multidimensional distributions the results of simulation seems to indicate quite a different behavior.</p>
<p>The simulations are done for dimensions (column &quot;dim&quot; in the table) 2,5,10,15,20,50,100 and samples number 1000,10 000, 100 000, 1 000 000 (column &quot;sample size&quot;) and each simulation was repeated 100 times (i.e. all simulation answers are averages over 100 samples).</p>
<p>When dimension grows the answer tends to 0.25 as proposed by Neal Young for sphere and torus , but for the  uniform and Gaussian it is not the case,  moreover the stabilization of the answers with sample size does not seen to happen for dimensions 50 and above.  The column &quot;Theory&quot; is answer proposed by Neal Young for corresponding dimension. (For dimension 2 it is result described in David Eppstein answer).</p>
<p><a href="https://i.sstatic.net/FX2BA.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/FX2BA.png" alt="Simulation results" /></a></p>
<p>Notebooks with simulations can be found here:
<a href="https://www.kaggle.com/alexandervc/study-connected-components-of-random-nn-graph?scriptVersionId=48963913" rel="nofollow noreferrer">https://www.kaggle.com/alexandervc/study-connected-components-of-random-nn-graph?scriptVersionId=48963913</a>
(V20 - torus V18,19 - sphere, V16 - uniform, V14,15 - Gaussian).
Simulation is done using GPU. Each notebook saves statistics of simulation to csv files - which are available in &quot;output&quot; section.
Summary table above can be found at the &quot;Data&quot; section of kaggle dataset environment:
<a href="https://www.kaggle.com/alexandervc/hubness-for-high-dimensional-datasets?select=NN_graphs_connected_components.xlsx" rel="nofollow noreferrer">https://www.kaggle.com/alexandervc/hubness-for-high-dimensional-datasets?select=NN_graphs_connected_components.xlsx</a>
Check that GPU and CPU libraries produce the same distribution is given e.g. here:
<a href="https://www.kaggle.com/alexandervc/knn-graphs-study-degree-distribution?scriptVersionId=48227434" rel="nofollow noreferrer">https://www.kaggle.com/alexandervc/knn-graphs-study-degree-distribution?scriptVersionId=48227434</a> (though due to precision issues the graphs themselves can be slightly different <a href="https://www.kaggle.com/alexandervc/benchmark-knn-graph-gpu-rapids-vs-sklearn?scriptVersionId=47940946" rel="nofollow noreferrer">https://www.kaggle.com/alexandervc/benchmark-knn-graph-gpu-rapids-vs-sklearn?scriptVersionId=47940946</a> ).</p>
<p>==================================================</p>
<p><strong>Older version</strong> (older version is mostly irrelevant - since it was not recognized that discrepancy appears for uniform and Gauss, but for sphere and torus it is Okay. But may be useful in respect that several checks that different implementation of K-NN graph construction seems to produce same results).</p>
<p>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 discrepancy appears for higher dimensions.</p>
<p>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.</p>
<p>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.</p>
<p>So let me share some more simulations results, self-checks, comments and scripts.</p>
<p>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...</p>
<p>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.</p>
<p>Dimension = 10 , Theory Percent 26.067</p>
<pre><code>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
</code></pre>
<p>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:
<a href="https://www.kaggle.com/alexandervc/connected-components-knn-graph-v010-rapids-knn?scriptVersionId=38115858" rel="nofollow noreferrer">https://www.kaggle.com/alexandervc/connected-components-knn-graph-v010-rapids-knn?scriptVersionId=38115858</a>
The simulation is using GPU package RAPIDs based on Facebooks FAISS <a href="https://engineering.fb.com/data-infrastructure/faiss-a-library-for-efficient-similarity-search/" rel="nofollow noreferrer">https://engineering.fb.com/data-infrastructure/faiss-a-library-for-efficient-similarity-search/</a>
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.</p>
<p><strong>What is subtle here:</strong>
It is known that GPU is single precision, while CPU is double precesion - and surprisingly enough that causes small difference in graphs produced. <a href="https://github.com/rapidsai/cuml/blob/34d473669bf323c26a38eff9e9a2c0454827688a/notebooks/nearest_neighbors_demo.ipynb" rel="nofollow noreferrer">(It is known)</a>.
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.</p>
<p>Dimension = 50 , Theory Percent about 25</p>
<pre><code>Sample Size   Percent by Simulation
1 000         16.694
10 000        15.6265
100 000       15.05882
1 000 000     14.834492
</code></pre>
<p>Notebook: <a href="https://www.kaggle.com/alexandervc/connected-components-knn-graph-v010-rapids-knn?scriptVersionId=38115858" rel="nofollow noreferrer">https://www.kaggle.com/alexandervc/connected-components-knn-graph-v010-rapids-knn?scriptVersionId=38115858</a></p>
<p>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.</p>
<p><strong>What is subtle here:</strong> See above</p>
<p>Dimension = 20 , Theory Percent about 25.19</p>
<pre><code>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
</code></pre>
<p>There is small increase with sample size, but theory is quite far...</p>
<p>Notebook up to 1 000 000 : <a href="https://www.kaggle.com/alexandervc/connected-components-knn-graph-v010-rapids-knn?scriptVersionId=37225738" rel="nofollow noreferrer">https://www.kaggle.com/alexandervc/connected-components-knn-graph-v010-rapids-knn?scriptVersionId=37225738</a>
Notebook  for 10 000 000 : <a href="https://www.kaggle.com/alexandervc/connected-components-knn-graph-v010-rapids-knn?scriptVersionId=37148875" rel="nofollow noreferrer">https://www.kaggle.com/alexandervc/connected-components-knn-graph-v010-rapids-knn?scriptVersionId=37148875</a></p>
<p>Dimension = 5 (100 times averaged)
Percent Theory = 27.8867</p>
<pre><code>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 
</code></pre>
<p>1e7 - time: 446.144 seconds - per 1 run
1e6 - time: 26.1098 seconds  - per 1 run</p>
<p><strong>What is subtle here:</strong> 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.</p>
<p>Here is notebook with speed comparaison:
<a href="https://www.kaggle.com/alexandervc/compare-nn-graph-speed-sklearn-vs-gpu-rapids" rel="nofollow noreferrer">https://www.kaggle.com/alexandervc/compare-nn-graph-speed-sklearn-vs-gpu-rapids</a></p>
<p>PS</p>
<p>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.</p>

CC BY-SA 4.0
2020-07-11T20:19:56.630Z
@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)  -  https://www.kaggle.com/alexandervc/rapids-knn-cugraph-test-9b97f4?scriptVersionId=38511131 - please look tables at the end of notebook ,  I will report more later ...
Alexander Chervov
CC BY-SA 4.0
2020-07-12T07:55:06.627Z
@NealYoung Sphere Dimension 50 , gives coefficient 25.05 - which finely agrees, with you prediction that results should be above 25 for any dimension. https://www.kaggle.com/alexandervc/rapids-knn-cugraph-test-9b97f4?scriptVersionId=38555171
Alexander Chervov
CC BY-SA 4.0
2020-07-13T02:29:54.497Z
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.
Neal Young
CC BY-SA 4.0
2020-07-13T08:09:45.407Z
@NealYoung sphere dimensions 500 and 5000 also give coefficent around 25 within std accuracy. For 500  https://www.kaggle.com/alexandervc/rapids-knn-cugraph-test-9b97f4?scriptVersionId=38615907 and for 5000 https://www.kaggle.com/alexandervc/rapids-knn-cugraph-test-9b97f4?scriptVersionId=38621415
Alexander Chervov
CC BY-SA 4.0
2020-07-14T17:30:12.703Z
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. https://www.kaggle.com/alexandervc/rapids-knn-cugraph-test-9b97f4?scriptVersionId=38673908
Alexander Chervov
CC BY-SA 4.0
2020-12-10T15:44:50.657Z
@NealYoung just updated simulation result - sphere and torus cases quite confirm your proposal
Alexander Chervov
CC BY-SA 4.0
2020-07-07T17:34:15.310Z
2020-12-10T15:55:45.017Z
2020-12-10T15:55:45.017Z
Alexander Chervov
<p>For <span class="math-container">$n$</span> uniformly random points in a unit square the number of components is</p>

<p><span class="math-container">$$\frac{3\pi}{8\pi+3\sqrt{3}}n+o(n)$$</span></p>

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

<p>For points in any fixed higher dimension it is <span class="math-container">$\Theta(n)$</span>; I don't know the exact constant of proportionality but the paper describes how to calculate it.</p>

CC BY-SA 4.0
2020-06-13T10:29:13.807Z
Thank you very much ! Your value is  0.3107524484437158 , so simulation above quite perfectly corresponds to it . What is Theta(n) ?
Alexander Chervov
CC BY-SA 4.0
2020-06-13T16:32:45.470Z
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.)
David Eppstein
CC BY-SA 4.0
2020-06-13T18:53:32.853Z
At the end of the paper you write: "It would be of interest to estimate some further properties of the NNG, such as the distribution of component sizes and the expected size of the largest component ". Are you aware of any result in that direction and related ? I mean since paper is 1997, it might be something more is known in 2020.
Alexander Chervov
CC BY-SA 4.0
2020-06-14T04:25:56.027Z
Hi David, do you know if the natural generalization to k nearest neighbor graphs have been looked at?
Sandeep Silwal
CC BY-SA 4.0
2020-06-14T05:20:44.513Z
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, https://doi.org/10.1007/s00453-007-9040-7
David Eppstein
CC BY-SA 4.0
2020-06-14T23:29:17.943Z
I see, thank you. What about the following question: In the same setting as above where we pick n points uniformly from the unit square, for what value of k is there only one (or O(1) ) connected component? Must k be Omega(n) ?
Sandeep Silwal
CC BY-SA 4.0
2020-06-15T06:59:53.920Z
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)$.
David Eppstein
CC BY-SA 4.0
2020-06-18T19:43:57.937Z
@DavidEppstein In   "k-Nearest-Neighbor Clustering and Percolation Theory"
they suggest that " k_2 = 3 , do not you know has this been proved  or discussed somewhere ?
Alexander Chervov
CC BY-SA 4.0
2020-06-18T19:57:34.063Z
I found that: Percolation in the k-nearest neighbor graph
P Balister, B Bollobás - Manuscript. Available from http …, 2008 . But that still 2008
Alexander Chervov
CC BY-SA 4.0
2020-06-18T19:57:43.093Z
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.
David Eppstein
CC BY-SA 4.0
2020-06-18T19:58:38.817Z
Thank you very much !
Alexander Chervov
CC BY-SA 4.0
2020-06-21T20:59:34.307Z
Any idea about the second order term ? To explain deviation between Neal Young and simulation it should be decreasing very very slowly with respect to "N" may be  (-1/Nlog(N)) at least for high DIMENSION
Alexander Chervov
CC BY-SA 4.0
2020-06-21T23:55:53.843Z
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.
David Eppstein
CC BY-SA 4.0
2020-06-13T04:33:23.697Z
2020-06-13T04:33:23.697Z
David Eppstein
<p>Let us sample some big number N points randomly uniformly on <span class="math-container">$[0,1]^d$</span>. 
Consider  <a href="https://en.wikipedia.org/wiki/Nearest_neighbor_graph" rel="noreferrer">1-nearest neighbor graph</a> based on such data cloud. (Let us look on it as UNdirected graph).</p>

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

<p>Simulation below suggest 31% for d=2,   20% for d=20, etc:</p>

<pre><code>Percent Dimension:
31 2
28 5
25 10
20 20
15 50
13 100
10 1000
</code></pre>

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

<p>If someone can comment on more general questions here: 
<a href="https://mathoverflow.net/q/362721/10446">https://mathoverflow.net/q/362721/10446</a>
that would be greatly appreciated.</p>

<pre><code>!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)
</code></pre>

CC BY-SA 4.0
2020-06-14T15:09:33.987Z
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.
Neal Young
<graph-theory><random-graphs>
Number of connected components of a random nearest neighbor graph?
Alexander Chervov
# Number of connected components of a random nearest neighbor graph?

(Asked by: username_0 on 2020-06-12)

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)

**Tags:** graph-theory, random-graphs

> **Comments:**
> **username_1**: 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.  

---

## 3 answers

### Answer by username_3 on 2020-06-13. Score: 17 | Accepted answer

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.

> **Comments:**
> **username_0**: Thank you very much ! Your value is  0.3107524484437158 , so simulation above quite perfectly corresponds to it . What is Theta(n) ?  
> **username_3**: 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.)  
> **username_0**: At the end of the paper you write: "It would be of interest to estimate some further properties of the NNG, such as the distribution of component sizes and the expected size of the largest component ". Are you aware of any result in that direction and related ? I mean since paper is 1997, it might be something more is known in 2020.  
> **username_2**: Hi David, do you know if the natural generalization to k nearest neighbor graphs have been looked at?  
> **username_3**: 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, https://doi.org/10.1007/s00453-007-9040-7  
> **username_2**: I see, thank you. What about the following question: In the same setting as above where we pick n points uniformly from the unit square, for what value of k is there only one (or O(1) ) connected component? Must k be Omega(n) ?  
> **username_3**: 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)$.  
> **username_0**: @DavidEppstein In   "k-Nearest-Neighbor Clustering and Percolation Theory" they suggest that " k_2 = 3 , do not you know has this been proved  or discussed somewhere ?  
> **username_0**: I found that: Percolation in the k-nearest neighbor graph P Balister, B Bollobás - Manuscript. Available from http …, 2008 . But that still 2008  
> **username_3**: 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.  
> **username_0**: Thank you very much !  
> **username_0**: Any idea about the second order term ? To explain deviation between username_1 and simulation it should be decreasing very very slowly with respect to "N" may be  (-1/Nlog(N)) at least for high DIMENSION  
> **username_3**: 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.  

---

### Answer by username_1 on 2020-06-13. Score: 11

*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](https://en.wikipedia.org/wiki/Volume_of_an_n-ball#The_volume) 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](https://en.wikipedia.org/wiki/Gamma_function#Stirling%27s_formula)).

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](https://math.stackexchange.com/a/162873/48535) 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](https://www.wolframalpha.com/input/?i=1+%2F%284-+2+BetaRegularized%5B3%2F4%2C+%28k%2B1%29%2F2%2C+1%2F2%5D%29+for+k%3D1%2C2%2C3%2C4%2C5%2C6%2C7%2C8%2C9%2C10%2C11%2C12%2C13%2C14%2C15%2C16%2C17%2C18%2C19%2C20).

---

$$
\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?

> **Comments:**
> **username_2**: very insightful answer, thank you :)  
> **username_0**: Thank you very much ! I will double check my simulation  
> **username_1**: @AlexanderChervov: If the analysis above is right, the number has to be between 1/4 and 1/2 for every dimension because $\mu_k$ is the volume of the union of two intersecting balls of volume $\beta_k$ each (so $\beta_k/(2\mu_k)$ has to be between 1/4 and 1/2).  So either the analysis or your code is wrong.  Re: the code, are you sure it's using the $\ell_2$ metric?  Could precision issues somehow be a problem?  One sanity check would be to check that for every vertex $u$ in the directed graph, the path $u\rightarrow v \rightarrow w$ (of nearest neighbors) has $d(v, w) \le d(v, u)$.  
> **username_1**: Also I'm fairly sure that $\lim_{k\rightarrow\infty} \beta_k/(2\mu_k) = 1/4$, because the  the volume of the "cap" becomes insignificant relative to the full ball as the dimension grows, so $\mu_k \sim 2\beta_k$.  
> **username_1**: FYI I sanity-checked your code and didn't find any mistake or obvious inconsistency.  
> **username_0**: Thank you for your remarks. I made preliminary checks, changed igraph to snap , got the same result . A friend of mine, using GPU studied greater size, like 1e6-7 points, he sees some increase of the answers for dim = 20 , but it is very low. Latter we will check a bit more.  
> **username_0**: # Dim = 20   # we are quite far from theory - 25.1904  
> **username_0**: n: 10**4 points -  mean percent : 20.486, time: 19.92 s  
> **username_0**: n: 10**5, mean percent: 20.897, time: 60.168 s  
> **username_0**: n: 2*10**5 points , mean percent: 21.044, time: 148.198 s  
> **username_0**: 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.  
> **username_0**: Later I will report on further simulations  
> **username_0**: Simulated dimension 50, it seems  percent is NOT growing with respect to sample size. n=100K 15.059 n=1Million 14.834 , with 100 times simulated and small std. Later I will report more. Need to check certain things  
> **username_1**: 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.  
> **username_0**: Thank you for your comment ! I have posted in answer my simulations with links and some comments (just in case). The main problem as it seems to me - for high dimensions - see D = 50 ( second example ). I will try to think more on that ....  
> **username_1**: @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](https://www.wolframalpha.com/input/?i=N%5B%288%2FPi%29%5E%28k%2F2%29+Gamma%28k%2F2%2B1%29%5D+where+k%3D14).  But perhaps refining the analysis could give tighter bounds.  
> **username_1**: @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.)  
> **username_1**: 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}$.  
> **username_0**: @NealYoung so how large should be "n" ? about ?  
> **username_0**: Python sklearn built-in functions should support several metrics: https://scikit-learn.org/stable/modules/generated/sklearn.neighbors.DistanceMetric.html?highlight=distancemetric%20class# But it seems you consider the torus (periodic) not the uniform cube , let me think how that can be done  
> **username_1**: Let us [continue this discussion in chat](https://chat.stackexchange.com/rooms/110359/discussion-between-neal-young-and-alexander-chervov).  
> **username_0**: What I can do just use take points in (S^1 \in R^2)^dimension - that corresponds to yours metric, is it correct ?  
> **username_0**: If you have some time please take a look: https://cstheory.stackexchange.com/questions/47957/power-law-for-degree-distribution-of-random-knn-graphs  
> **username_1**: @AlexanderChervov I took a look and I don't know.. My guess would be that the maximum in-degree of any node is bounded by some function f(d, k) that depends only on d and k, independently of the number of points or how they are located (randomly, adversarially, ...).  Can you construct a counterexample to this guess?  Hopefully somebody who actually knows something about the literature (unlike me) will have something to say... :-)  
> **username_0**: @NealYoung thank you very much for your answer I will try to check.  
> **username_0**: @NealYoung I tried to check whether maximum is bounded - it seems it is not (at least for Gaussian, dim =50) - see https://stats.stackexchange.com/questions/500250/recognize-distribution-on-integers-0-mean-5-median-2-percentiles-75-5   The distribution seems to have some nice properties like std apprximately 2*mean, mean = percentile 75  

---

### Answer by username_0 on 2020-07-07. Score: 0

Update:
Not a answer but sharing simulation result confirming proposed answer by username_1 for certain cases.
username_1 in his answer above a proposed beautiful formula for the number of connected components of NN-graph in any dimension. Simulation results in the table below quite confirm his proposal for the case of points distributed uniformly on sphere and torus (probably any closed manifold). Though for points taken from the uniform (on cube) or Gaussian multidimensional distributions the results of simulation seems to indicate quite a different behavior.
The simulations are done for dimensions (column "dim" in the table) 2,5,10,15,20,50,100 and samples number 1000,10 000, 100 000, 1 000 000 (column "sample size") and each simulation was repeated 100 times (i.e. all simulation answers are averages over 100 samples).
When dimension grows the answer tends to 0.25 as proposed by username_1 for sphere and torus , but for the  uniform and Gaussian it is not the case,  moreover the stabilization of the answers with sample size does not seen to happen for dimensions 50 and above.  The column "Theory" is answer proposed by username_1 for corresponding dimension. (For dimension 2 it is result described in username_3 answer).

Notebooks with simulations can be found here:
https://www.kaggle.com/alexandervc/study-connected-components-of-random-nn-graph?scriptVersionId=48963913
(V20 - torus V18,19 - sphere, V16 - uniform, V14,15 - Gaussian).
Simulation is done using GPU. Each notebook saves statistics of simulation to csv files - which are available in "output" section.
Summary table above can be found at the "Data" section of kaggle dataset environment:
https://www.kaggle.com/alexandervc/hubness-for-high-dimensional-datasets?select=NN_graphs_connected_components.xlsx
Check that GPU and CPU libraries produce the same distribution is given e.g. here:
https://www.kaggle.com/alexandervc/knn-graphs-study-degree-distribution?scriptVersionId=48227434 (though due to precision issues the graphs themselves can be slightly different https://www.kaggle.com/alexandervc/benchmark-knn-graph-gpu-rapids-vs-sklearn?scriptVersionId=47940946 ).
==================================================
Older version (older version is mostly irrelevant - since it was not recognized that discrepancy appears for uniform and Gauss, but for sphere and torus it is Okay. But may be useful in respect that several checks that different implementation of K-NN graph construction seems to produce same results).
It is not an answer, but a comment to very interesting answer by username_1;
who proposed the beautiful formula for general dimensions; thus generalizing the beautiful formula  by username_3 and coauthors. The formula fits simulations very nicely for lower dimensions; however the discrepancy 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.

> **Comments:**
> **username_0**: @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)  -  https://www.kaggle.com/alexandervc/rapids-knn-cugraph-test-9b97f4?scriptVersionId=38511131 - please look tables at the end of notebook ,  I will report more later ...  
> **username_0**: @NealYoung Sphere Dimension 50 , gives coefficient 25.05 - which finely agrees, with you prediction that results should be above 25 for any dimension. https://www.kaggle.com/alexandervc/rapids-knn-cugraph-test-9b97f4?scriptVersionId=38555171  
> **username_1**: 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.  
> **username_0**: @NealYoung sphere dimensions 500 and 5000 also give coefficent around 25 within std accuracy. For 500  https://www.kaggle.com/alexandervc/rapids-knn-cugraph-test-9b97f4?scriptVersionId=38615907 and for 5000 https://www.kaggle.com/alexandervc/rapids-knn-cugraph-test-9b97f4?scriptVersionId=38621415  
> **username_0**: 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. https://www.kaggle.com/alexandervc/rapids-knn-cugraph-test-9b97f4?scriptVersionId=38673908  
> **username_0**: @NealYoung just updated simulation result - sphere and torus cases quite confirm your proposal  

---

