Integration of Booles and Gaussians in high dimensions

Question

Context

As a follow up of this question, I would like to predict the connectivity of the so-called cosmic web in arbitrary dimensions.

The connectivity $\kappa$ is defined as the number of ridges connecting a given maxima to its surrounding saddles, which in turn therefore corresponds to $2 n_{\rm max}/ n_{\rm saddle}$ where $n_{\rm max}$ and $n_{\rm saddle}$ are the total number of maxima and saddles in the field.

This is how the cosmic web looks like in 3D:

so the motivation is to predict how many filaments are connected to a given node from first principes.

Attempt

As a starting point I am restricting myself to a Gaussian random field. For such fields it can be shown that the following code computes this ratio for an arbitrary dimension $d$.

d=3; 
var = Table[ToExpression["x" <> ToString[i]], {i, d}];
mat = ConstantArray[1/(d (d + 2)) , {d, d}];
Do[mat[[i, i]] = 3/(d (d + 2)), {i, d}]; 
arg = Abs[Times @@ var] Abs[
    Product[var[[i]] - var[[j]], {i, d}, {j, i - 1}]];
pdf = MultinormalDistribution[Table[0, {Length[mat]}], mat];
fun1 = Function[var // Evaluate, 
   Times @@ (Boole[# < 0] & /@ var) arg // Evaluate];
fun2 = Function[var // Evaluate, 
   Boole[Last[var] > 0] Times @@ (Most@(Boole[# < 0] & /@ var)) arg //
     Evaluate];

2   Integrate[fun2 @@ var PDF[pdf, var],
   Sequence @@ Map[{#, -Infinity, Infinity} &, var]]/
  Integrate[fun1 @@ var PDF[pdf, var],
   Sequence @@ Map[{#, -Infinity, Infinity} &, var]]

For instance it returns 4 for $d=2$ and

2 (18 Sqrt[2] + 29 Sqrt[3]) / (-18 Sqrt[2] + 29 Sqrt[3]) = 6.11

for $d=3$. It is strikingly close to a cubic face centred lattice, but with some level of impurity to the crystal.

for $d>3$ I have resorted to NIntegrate

2 d  NIntegrate[fun2 @@ var PDF[pdf, var] // Evaluate,
   Sequence @@ Map[{#, -Infinity, -1, 0, 1, Infinity} &, var] // 
    Evaluate, PrecisionGoal -> 20, MinRecursion -> 8]/
  NIntegrate[fun1 @@ var PDF[pdf, var] // Evaluate,
   Sequence @@ Map[{#, -Infinity, -1, 0, 1, Infinity} &, var] // 
    Evaluate,
   PrecisionGoal -> 20, MinRecursion -> 8]

and I find 8.35 10.75 13.23 for $d$=4,5,6

Question

Is is possible to compute the ratio analytically for $d>3$?

FYI, for $d=4$ the integrant of the numerator looks like this

      int = fun2 @@ var PDF[pdf, var]; int

Alternative Question

What options should one use with NIntegrate so that the code works for $d>6$.

Scientifically, it would be of interest to go up to d=11 in the context of so-called landscape inflation.

As a test, I have tried the following brute force Monte Carlo:

 d=5;
 dat2 = ParallelTable[dat = RandomVariate[pdf, 150000];
 {2 d (fun2 @@ # & /@ dat // Total), (fun1 @@ # & /@ dat // 
   Total)}, {25}]; // AbsoluteTiming
 Mean@ dat2 // First[#]/Last[#] &

which is embarrassingly parallel but not very accurate.

So far I have

$\kappa^d$ = $4$, $6.11$, $8.35$, $10.73$, $13.23$, $15.85$ for $d=$ $2$, $3$, $4$, $5$, $6$, $7$ resp.

From inspection of $d\le 7$ we conjecture that it is closely approximated by

$$ \kappa^d= 2d+\left(\frac{2d-4}{7}\right)^{7/4}. $$

PS: I have now also asked this question on math.stackexchange

Answer 1

Some partial results.

So I've been working on this problem in my spare time, and I wanted to post some partial results I got should it help someone else (and, of course, OP) to work through the problem.

First of all, a more-or-less compact expression for the integral OP is after:

-2 d Integrate[
 Times @@ Array[λ, d] Abs[Product[λ[j] - λ[i], {i, 1, d}, {j, i + 1, d}]] Exp[-1/2 Sum[λ[i] λ[j] (-d/2), {i, 1, d}, {j, 1, d}] - 1/2 Sum[λ[i]^2 1/2 d (2 + d), {i, 1, d}]]
     , Evaluate[Sequence @@ Table[{λ[i], -∞, 0}, {i, 1, d - 1}]], {λ[d], 0, ∞}]/
Integrate[
 Times @@ Array[λ, d] Abs[Product[λ[j] - λ[i], {i, 1, d}, {j, i + 1, d}]] Exp[-1/2 Sum[λ[i] λ[j] (-d/2), {i, 1, d}, {j, 1, d}] - 1/2 Sum[λ[i]^2 1/2 d (2 + d), {i, 1, d}]]
     , Evaluate[Sequence @@ Table[{λ[i], -∞, 0}, {i, 1, d}]]]

Integrating numerically using the options Method -> "GaussKronrodRule", PrecisionGoal -> 20, AccuracyGoal -> 20, WorkingPrecision -> 30, MaxRecursion -> 50, I found $$ \kappa^{(4)}=\color{red}{8.34922775472958914088481}88794 $$ of which I am pretty sure the red figures are significative. Perhaps a few more figures are as well.

It would be great to evaluate this integral to a better precision, and use RootApproximant to guess the analytic result. If, as in the $d=2,3$ cases, the result is algebraic, RootApproximant will be able to reconstruct it from the numerical value if given enough correct digits. My attempts at doing this appear to have failed, so my guess is that $\kappa^{(d)}$ is not algebraic for $d\ge4$. Oh well.

If it turns out that the analytic form of $\kappa^{(4)}$ can be obtained, and so can the next few values, perhaps FindSequenceFunction will be able to guess the value of $\kappa^{(d)}$ for arbitrary $d$. Wouldn't it be nice?

A possible approach to improve the convergence is to enforce some change of variables. For now, I managed to prove that the change of variables $\vec\lambda\to A\vec\lambda$ diagonalises the exponent, where $A$ is the (orthogonal) matrix given by

A = Join[Table[PadRight[Join[ConstantArray[-(1/Sqrt[j (1 + j)]), j], {Sqrt[j/(1 + j)]}], d], {j, 1, d - 1}], {ConstantArray[1/Sqrt[d], d]}];

In particular, the inverse of the matrix $M$ in the OP satisfies

ConstantArray[-d/2, {d, d}] + 1/2 d (2 + d) IdentityMatrix[d] == Transpose[A].DiagonalMatrix[Join[ConstantArray[1/2 d (2 + d), d - 1], {d}]].A

so that the eigenvalues are $d-1$ copies of 1/2 d (2 + d), and one copy of d.

My hope is that this basis allows us to factorise the integral and evaluate it analytically for arbitrary $d$. I haven't had the time to explore this yet.

Finally, let me mention that I also evaluated $\kappa^{(d)}$ numerically using the options Method -> "MonteCarlo" (with increasing number of points to estimate the error), and I found $$ \begin{aligned} \kappa^{(8)}&=18.7(1)\\ \kappa^{(9)}&=21.4(2)\\ \kappa^{(10)}&=24.4(2)\\ \kappa^{(11)}&=27.4(2) \end{aligned} $$

The evaluation of $\kappa^{(d)}$ for higher $d$ is straightforward: just use the integral expression I wrote above, but change Integrate to NIntegrate and use the option Method -> "MonteCarlo". You will need to increase the number of points though: for $d=11$ I had to use MaxPoints -> 700000000. For higher $d$, you'll probably need to use some more. And to estimate the error, you can evaluate NIntegrate a few times, and take the mean and standard deviation. At this point, the evaluation of $\kappa^{(d)}$ is just a matter of patience. I guess at some point it will be necessary to increase the working precision too. I leave this to you, if you feel more values of $\kappa^{(d)}$ are necessary.

Finally, let me mention that the results so far are all in fact consistent with OP's conjecture:

Show[
     ListPlot[Transpose[{Range[2, 11], {4, 6.11, 8.35, 10.73, 13.23, 15.85, 18.7, 21.4, 24.4, 27.4}}]],
     Plot[{2 d, 2 d + ((7 d - 14)/25)^(7/4)}, {d, 2, 11}]
]

I'll update this post as I get more results.

Answer 2

Through some indirect theoretical consideration it turns out it can be shown that the connectivity for $d=4$ is in fact symbolic and equal to

 -((200 \[Pi])/(114-75 \[Pi]+100 cot^-1(2))) = 8.3492277547295891408846033687519774106386215942178

So I can confirm that @accidentalFourierTransform's result is correct up to 22 digits :-)

The sketch of the proof is given here