How to speed up this high-dimensional definite integral?

Question

Examples

➤ 1.

$$I=\iiint\limits_D \dfrac{{\rm d}x\ \!{\rm d}y\ \!{\rm d}z}{(x+y+z)^p}~,$$

where $$D=\{~(x,y,z)~\colon~x>0~,y>0~,z>0~,x+y+z<1~\}~.$$

MMA gives: $~I=-\dfrac1{2(p-3)}~.$

➤ 2.

$$I=\idotsint\limits_{~~V} \dfrac{{\rm d}x_1\ \!{\rm d}x_2\ \!\cdots\ \!{\rm d}x_n}{(x_1+x_2+...+x_n)^p}~,$$

where $$V=\{~(x_1,x_2~,\dots,x_n)~\colon~x_1>0~,x_2>0~~,\dots,x_n>0~,x_1+x_2+...+x_n<1~\}~.$$

The integral converges when $p<n$.

Verified the case of $n=1～5$, conjecture by induction that $$I=\frac{1}{(n-1)!~(n-p)}~.$$

---

Conclusion

Special thanks to @Pillsy's formula and @CarlWoll's simplification!

In:

n = 20;
Integrate[y^(n-1-p), {y,0,1}, x∈Simplex[n-1]] [[1]] // AbsoluteTiming

Out:

$\left\{0.363959,\frac1{121645100408832000 (20-p)}\right\}$

---

NOTE: The following content is no longer worth browsing.

---

Code

In:

int[var_] :=
    Block[
        var,
        fun = Power[Plus @@ var, p]~Power~-1;

        con0 = <|1 -> (StringRiffle[# > 0 & /@ var, "&&"] // ToExpression), 
                 2 -> Plus @@ var < 1|>;

        con = StringRiffle[Values@con0, "&&"] // ToExpression;

        Column[{fun, con}] // TraditionalForm // Print;

        (result = Integrate[fun, var ∈ ImplicitRegion[con, Evaluate@var], 
                      Assumptions -> p < Length@var] // Together
        ) // TraditionalForm
    ]

int@{x, y, z}
int@{x, y, z, u, v, w} // Timing

Out:

$(x+y+z)^{-p} \\ x>0\land y>0\land z>0\land x+y+z<1 \\ -\dfrac1{2(p-3)}$

$(u+v+w+x+y+z)^{-p} \\ x>0\land y>0\land z>0\land u>0\land v>0\land w>0\land u+v+w+x+y+z<1 \\ \left\{29.6719,-\dfrac1{120(p-6)}\right\}$

---

Timing

$\begin{array}{rr} n & t(s) \\ \hline 1 & 0.328 \\ 2 & 0.656 \\ 3 & 1.391 \\ 4 & 3.297 \\ 5 & 16.828 \\ 6 & 29.641 \\ 7 & 69.766 \end{array}$

Tests on:
Mathematica 11.1.1.0 (←2017-05-16),
Windows 64-bit 10.0.10586 (←2015-11-12),
AMD A8-5600K (←2012-09-26)

---

Aiming

Code of fast and efficient algorithm.

Wish: regardless of how large is $n$, the computation time is always less than 10s.

---

I would be thankful, if you could improve or rewrite any part of the code.

I would be admired, if you could write a library/macro that can output very fast and can apply to many situations, whereas even do not depend on the original Integrate function.

Question END.

---

---

The following has nothing to do with the post topic, just recording my experiences of fixing this —

Your post appears to contain code that is not properly formatted as code.

In my case, all the problems stem from LaTeX MathJax (comment from @anderstood).
① Custom commands are not allowed to submit (though they are rendered normally).
② Need to delete ⏎ after \\ in some places.

Answer 1

I was able to speed things up incredibly by observing that the integrand only depends on the sum $ x_1 + x_2 + \ldots + x_n $, and we can take advantage of this by introducing a gratuitous integration over a new variable and throwing in a Dirac $ \delta $ "function", which satisfies the following relation:

$$ \int_a^b dx\,f(x)\,\delta(x-c) = \begin{cases} f(c) & a < c < b \\ 0 & c < a \text{ or } b < c \end{cases} $$

With a little coaxing, you can reproduce this in Mathematica, using DiracDelta:

Integrate[f[x]*DiracDelta[x - c], {x, a, b},
   Assumptions -> {a < b, #}] & /@ {
   a < c < b,
   c < a || b < a
  }
(* {f[c], 0} *)

Now, use a specific application of the defining relation. For $ 0 < x_1 + \ldots + x_n < 1 $, we have

$$ \frac{1}{\left(x_1 + \ldots + x_n\right)^{\,p}} = \int_0^1 dy \frac{\delta\left(y - (x_1 + \ldots + x_n\right)}{y^p} $$

We can insert this into the original expression, and swap around the order of the integration.

$$ \idotsint\limits_{V} \frac{dx_1\cdots\,dx_n}{(x_1 + \ldots + x_n)^p } = \int_0^1 \frac{dy}{y^p} \int_0^\infty dx_1 \dots \int_0^\infty dx_n \,\,\, \delta\left(y - (x_1 + \ldots + x_n\right) $$

Now we can do the integral over the $x_i$s separately

deltaFunctionIntegral[n_Integer?Positive, y_] :=
 With[{vars = Array[\[FormalX], n]},
  Apply[
   Integrate[DiracDelta[y - Total[vars]], ##, 
     Assumptions -> Positive[y]] &,
   Map[{#, 0, Infinity} &, vars]]];

Mathematica will blaze through it:

Table[deltaFunctionIntegral[n, y] // AbsoluteTiming, {n, 1, 10}]
(* 
   {{0.024917, 1}, 
    {0.033459, y}, 
    {0.042226, y^2/2}, 
    {0.045033, y^3/6}, 
    {0.050282, y^4/24}, 
    {0.05671, y^5/120}, 
    {0.063595, y^6/720}, 
    {0.073157, y^7/5040}, 
    {0.080864, y^8/40320}, 
    {0.086789, y^9/362880}}
*)

Now we can do the integral with this factor and things remain speedy:

Table[
  With[{factor = deltaFunctionIntegral[n, y]},
    Integrate[factor/y^p, {y, 0, 1}]] // AbsoluteTiming
  {n, 1, 10}]
(* 
   {{0.288039, ConditionalExpression[1/(1 - p), Re[p] < 1]}, 
    {0.300244, ConditionalExpression[1/(2 - p), Re[p] < 2]}, 
    {0.314402, ConditionalExpression[1/(2 (3 - p)), Re[p] < 3]}, 
    {0.32506, ConditionalExpression[1/(6 (4 - p)), Re[p] < 4]},
    {0.318678, ConditionalExpression[1/(24 (5 - p)), Re[p] < 5]}, 
    {0.339789, ConditionalExpression[1/(120 (6 - p)), Re[p] < 6]}, 
    {0.346237, ConditionalExpression[1/(720 (7 - p)), Re[p] < 7]}, 
    {0.35325, ConditionalExpression[1/(5040 (8 - p)), Re[p] < 8]}, 
    {0.365873, ConditionalExpression[1/(40320 (9 - p)), Re[p] < 9]}, 
    {0.367103, ConditionalExpression[1/(362880 (10 - p)), Re[p] < 10]}}
*)

I know how to show that deltaFunctionIntegral[n, y] gives $ y^{n-1}/(n-1)! $ outside of Mathematica for general $ n $, but I don't know how to coax the result out of Mathematica yet.

Answer 2

This is what I get in my 7-years old Intel I3 in my macOS 10.13.2 and MMA ver 11.2:

Table[AbsoluteTiming[vars = Array[x, n]; 
Integrate[(Apply[Plus, vars])^-p, vars \[Element] 
ImplicitRegion[And @@ (# > 0 & /@ vars) && Apply[Plus, vars] < 1,Evaluate@vars]]], {n, 7}]

$\left( \begin{array}{cc} 0.488669 & \text{ConditionalExpression}\left[\frac{1}{1-p},\Re(p)<1\right] \\ 0.978526 & \text{ConditionalExpression}\left[\frac{\frac{1}{2-p}-1}{p-1},\Re(p)<2\right] \\ 2.07194 & \text{ConditionalExpression}\left[\frac{\frac{1}{3-p}-\frac{p}{2}}{p^2-3 p+2},\Re(p)<3\right] \\ 4.77045 & \text{ConditionalExpression}\left[\frac{p^2-2 p-\frac{6}{4-p}+3}{-6 p^3+36 p^2-66 p+36},\Re(p)<4\right] \\ 20.0045 & \text{ConditionalExpression}\left[-\frac{1}{24 (p-5)},\Re(p)<5\right] \\ 33.4052 & \text{ConditionalExpression}\left[-\frac{1}{120 (p-6)},\Re(p)<6\right] \\ 84.0427 & \text{ConditionalExpression}\left[\frac{\frac{720}{7-p}-p ((p-7) p ((p-7) p+28)+252)}{720 (p-6) (p-5) (p-4) (p-3) (p-2) (p-1)},\Re(p)<7\right] \\ \end{array} \right)$

Answer 3

(This answer is just a simplification of @Pillsy's answer)

You could make the change of variables $x_i \to y x_i$ to convert the integral to $\underset{x\in \Delta^{n-1}}{\int } y^{n-1} dx$ (here $\Delta^{n-1}$ is the usual nomenclature for Simplex[n-1]) which evaluates very quickly. For example:

Integrate[y^99, x ∈ Simplex[99]] //AbsoluteTiming

{0.02938, y^99/ 933262154439441526816992388562667004907159682643816214685929638952175999932299156089414639761565182862536979208272237582511852109168640000000000000000000000}

By comparison, deltaFunctionIntegral takes about 1000 times longer:

deltaFunctionIntegral[100, y] //AbsoluteTiming

{31.0896, y^99/ 933262154439441526816992388562667004907159682643816214685929638952175999932299156089414639761565182862536979208272237582511852109168640000000000000000000000}

Answer 4

Old school, directly construct the integration limits instead of using a region proves to be just slightly faster.

n = 5
var = Symbol["a" <> ToString[#]] & /@ Range[n]
fun = 1/(Total@var)^p
Integrate[fun, Evaluate@Sequence @@ MapThread[{#1, 0, 1 - #2} &,
     {var, Accumulate[Prepend[var[[;; -2]], 0]]}], 
  Assumptions -> p < Length@var] // Timing

{9.90606, 1/(120 - 24 p)}

n=6

{29.5466, 1/(720 - 120 p)}

It doesn't seem to make much difference in timing if you leave out the assumption on p and let it generate a conditional.