How to speed up this high-dimensional definite integral?
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)$
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.
(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}