2
$\begingroup$

I'm trying to understand large $d$ behavior for $F(p,d)$ with $1\le p\le 2$ and $d>10^6$

$$F(p, d)=\frac{1}{d}\frac{\sum_{ij}^d(i^{-p}-j^{-p})^2}{\left(\sum_i^d i^{-p}\right)^2}$$

The following code works, but doesn't scale for large $d$, can someone see a way to get this to get this plotted for large $d$?

Clear[h, i, j];
var[p_, d_] := (
   h[i_] = i^-p;
   1/d NSum[(h[i] - h[j])^2, {i, 1, d}, {j, 1, d}]/
    NSum[h[i], {i, 1, d}]^2
   );
DiscretePlot[{var[1., d], var[1.5, d], var[2., d]}, {d, 1., 20, 1.}, 
 AxesLabel -> {"d", "value"}, PlotLegends -> {"p=1", "p=1.5", "p=2"}]

enter image description here

I'm interested in finding what the limit is, and how long it takes to get there as a function of $p$

Motivation: this gives variance of $f(x)=x_1/1^p + x_2/2^p + x_3/3^p+\ldots+x_d/d^p$ evaluated at a random point $x$ on a sphere of radius chosen such that expected value $E[f(x)]=1$ (variance formula from here). The case for $p=1$ is worked out here, I'm trying to get the $p>1$ case.

$\endgroup$

3 Answers 3

7
$\begingroup$

If you change NSum to Sum, then you get

var[p_, d_] := (h[i_] = i^-p;
   1/d Sum[(h[i] - h[j])^2, {i, 1, d}, {j, 1, d}]/Sum[h[i], {i, 1, d}]^2);

with the general result

var[p, d]
(* -((2 (HarmonicNumber[d, p]^2 - d HarmonicNumber[d, 2 p]))/(d HarmonicNumber[d, p]^2)) *)

The limit with $p=2$ as $d \rightarrow\infty$ is

Limit[-((2 (HarmonicNumber[d, 2]^2 - d HarmonicNumber[d, 4]))/(d HarmonicNumber[d, 2]^2)), d -> ∞]
(* 4/5 *)

Based on the comment from @UlrichNeuman a more general result can be obtained:

Limit[var[p, d], d -> ∞, Assumptions -> p > 1]
(* (2 Zeta[2 p])/Zeta[p]^2 *)

Using Mathematica 13.2.0.0 on Windows 10 I get the following for $0\leq p\leq 1$:

Limit[var[p, d], d -> ∞, Assumptions -> 0 <= p <= 1]

result for restricted range of p

I don't know why the rest of the range of $p$ is not included but the complete range can be covered with

Limit[var[p, d], d -> ∞, Assumptions -> 0 <= p <= 1/2]
(* 0 *)

and

Limit[var[p, d], d -> ∞, Assumptions -> 2/3 <= p <= 1]

Limit for 2/3 <= p <= 1

$\endgroup$
2
  • 2
    $\begingroup$ Perhaps Asymptotic[var[p, d], d -> Infinity, Assumptions -> 1 < p < 2] gives (2 Zeta[2 p])/Zeta[p]^2 might complete your answer $\endgroup$ Commented Jun 6, 2023 at 8:54
  • $\begingroup$ Thanks, @UlrichNeumann! $\endgroup$
    – JimB
    Commented Jun 6, 2023 at 17:17
4
$\begingroup$

If $F(p, d)=\frac{1}{d}\frac{\sum_{ij}^d(i^{-p}-j^{-p})^2}{\left(\sum_i^d i^{-p}\right)^2}$, then asymptotically $F(p,d) = \frac{2 \zeta (2 p)}{\zeta (p)^2} + O(1/d)$. So the asymptotic answer is $\lim_{d->\infty} F(p,d) = \frac{2 \zeta (2 p)}{\zeta (p)^2}$.

2 Zeta[2 p]/Zeta[p]^2 /. p -> 2

4/5

2 Zeta[2 p]/Zeta[p]^2 /. p -> 1.5

0.352277

$\endgroup$
3
$\begingroup$

This makes me think that one should consider using mathematics instead of Mathematica:

Compile[{{d, _Integer}, {p, _Real}},
   Block[{hh},
    hh = Range[1., d]^-p;
    1/d Total[Total[(hh - #)^2] & /@ hh]/Total[hh]^2
    ]][10^6, 1.5] // AbsoluteTiming

(*  {219.747, 0.352815}  *)

But maybe it's fast enough....

$\endgroup$
1
  • $\begingroup$ Interesting, plotting it for $d=10^5,p=2$ gives result suspiciously close to $4/5$ so there might be a simple closed form solution, posted on math.SE $\endgroup$ Commented Jun 6, 2023 at 1:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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