# Computing $\frac{1}{d}\int_1^d \int_1^d \mathrm{d}i\ \mathrm{d}j\ \left(f(i,s)-f(j,s)\right)^2$ for $d\to\infty$

I'm looking to compute the following integral

$$U(s)=\frac{1}{d}\int_1^d \int_1^d \mathrm{d}i\ \mathrm{d}j\ \left(f(i,s)-f(j,s)\right)^2$$

Where $$f(i,s)=i^{-p}\exp(-2s i^{-p}), s>1, p>1, d\approx \infty$$.

Code below takes 32 seconds and finds that for $$p=6/5$$, $$U(s)\approx \frac{5 \Gamma \left(\frac{7}{6}\right)}{12 \sqrt[3]{2} s^{7/6}}$$.

Is it possible to get a nice formula for generic $$p>1$$ in reasonable time?

ClearAll["Global*"];
$Assumptions = {d > 0, s > 0}; p = 1 + 1/5; h[i_] = (i + 1)^-p; f[i_, s_] = h[i] Exp[-2 s h[i]]; Asymptotic[Integrate[1/d (f[i, s] - f[j, s])^2, {i, 0, d}, {j, 0, d}], d -> Infinity]  Motivation: modeling of gradient descent, math.SO q ## 2 Answers Perhaps I misunderstood the question... Adding the constraint d>1 helps MMA to solve the integral h[i_] = (i )^-p; f[i_, s_] = h[i] Exp[-2 s h[i]]; int = Integrate[1/d (f[i, s] - f[j, s])^2, {i, 1, d}, {j, 1, d},Assumptions -> {s > 1, p > 1, d > 1}]  asy1 = Asymptotic[int, {d, Infinity, 1} ,Assumptions -> {s > 1, p > 1} ] // Simplify  Hope it helps! Adding this as an answer because it's kind of long for a comment, but it does not answer your question completely: First, is the lower bound supposed to be one or zero? I was able to evaluate a closed form solution for generic p >1 for at least just the integral w.r.t. j so at least one of the integrations can be done: ClearAll["Global*"];$Assumptions = {d > 0, s > 0, p > 1};
h[i_] = (i + 1)^-p;
f[i_, s_] = h[i] Exp[-2 s h[i]];
integrand = (f[i, s] - f[j, s])^2;
justJ = FullSimplify[Integrate[integrand, {j, 0, d}]]

(*(1/(16 p^3 s^2))(p (16 E^(-2 (1 + (1 + i)^-p) s) (1 + i)^-p p s -
16 (1 + d) E^(-2 ((1 + d)^-p + (1 + i)^-p) s) (1 + i)^-p p s +
16 d E^(-4 (1 + i)^-p s) (1 + i)^(-2 p) p^2 s^2 -
E^(-4 s) (-1 + p + 4 p s) + (1 + d) E^(-4 (1 + d)^-p s) (-1 +
p + 4 (1 + d)^-p p s)) -
16 E^(-2 (1 + i)^-p s) (1 + i)^-p p s ExpIntegralE[1 + 1/p,
2 s] + (-1 + p) ExpIntegralE[1 + 1/p, 4 s] +
16 (1 + d) E^(-2 (1 + i)^-p s) (1 + i)^-p p s ExpIntegralE[1 + 1/p,
2 (1 + d)^-p s] - (1 + d) (-1 + p) ExpIntegralE[1 + 1/p,
4 (1 + d)^-p s])*)



This took about 10s on my computer.

Edit: It takes only another ~15s to evaluate:

final = FullSimplify[Integrate[justJ, {i, 0, d}]]

(* (1/(p^2))(1 +
d)^(-2 p) (-2 d (1 + d)^(2 p) p ExpIntegralE[-1 + 1/p, 4 s] +
2 d (1 + d) p ExpIntegralE[-1 + 1/p, 4 (1 + d)^-p s] -
2 ((1 + d)^
p ExpIntegralE[1/p, 2 s] - (1 + d) ExpIntegralE[1/p,
2 (1 + d)^-p s])^2) *)



But I cannot get Limit[final/d, d -> Infinity] to evaluate. I am also not 100% sure if I've made any mathematical sin by doing the integrals first and then taking the limit as d->Infinity

Final Edit:

after going to get coffee, (so < 5 min), the limit is evaluated as:

Limit[final/d, d -> Infinity]

(*(2^(-3 + 2/p) s^(-2 + 1/p) (Gamma[2 - 1/p] - Gamma[2 - 1/p, 4 s]))/p*)



sorry for all of the iterative edits on this response.

• Adding a comment, if the lower integral bounds are supposed by 1 and not 0 (I used 0 in my response), then the code still runs but takes a total of about 10 minutes from start to finish to produce Limit[final/d, d -> Infinity] which ends up being (2^(-3 + 2/p) s^(-2 + 1/ p) (Gamma[2 - 1/p] - Gamma[2 - 1/p, 2^(2 - p) s]))/p
– ydd
Jun 7, 2023 at 18:14
• Are you sure the answers in comments are correct? I get nothing like them. Jun 7, 2023 at 18:33
• @1729taxi for the case where the lower bound is 0 yes, I just ran it again and checked the commented output and they are correct. I will check the lower bound = 1 case again now. What did you get?
– ydd
Jun 7, 2023 at 18:39
• A huge mess is what I get. No ExpIntegralE in either of the results you have in comments and the limit just goes off for ever. Jun 7, 2023 at 18:43
• You ran exactly the input I have in the answer, with the ClearAll["Global*"] and \$Assumptions` ? I just ran it again even after quitting and restarting mma and got the same output I have in the answer.
– ydd
Jun 7, 2023 at 18:50