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$

Question

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

Answer 1

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!

Answer 2

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.