# Plotting $\frac{1}{d}\frac{\sum_{ij}^d(i^{-p}-j^{-p})^2}{\left(\sum_i^d i^{-p}\right)^2}$ for large $d$

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"}]


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.

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]


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]


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

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

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....

• 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 Commented Jun 6, 2023 at 1:37