Numerical solution of high-dimensional integral involving positive-part function

Edit: Unfortunately I can not comment on your question. Perhaps someone can like my answer to get a reputation of 50. That's no joke!

Anyway, you may use a Gauss Laguerre quadrature or Gauss Hermite quadrature to calculate it. The quadrature rules are designed for integration kernels of following form:

Laguerre quadrature:

$$int_{0}^{infty} e^{-x} f(x),dx approx sum_{i=1}^n w_i f(x_i),$$

where

$$w_i = frac {x_i} {left(n + 1right)^2 left[L_{n+1}left(x_iright)right]^2}$$

and $L_n$ is the n-th Laguerre polynomial.

Hermite quadrature:

$$int_{-infty}^{+infty} e^{-x^2} f(x),dx approx sum_{i=1}^n w_i f(x_i),$$

where

$$w_i = frac {2^{n-1} n! sqrt{pi}} {n^2[H_{n-1}(x_i)]^2}$$

and $H_n$ is the n-th Hermite polynomial.

Note:

The given quadrature rules should rather be considered as an indication. They should easily extendable to multivariate functions.

You mentioned:

"I've been facing this integral for quite some time now, with the goal of disentangling the max-function and arriving at a way to perform the summations separately."

Perhaps you can elaborate more on your comment?

Let's consider the integral in 2D.

Note that the domain where $f(x_1,x_2) = x_1 + x_2 - K$ is positive lies to the right of the dashed line $x_1+x_2=K$ in the sketch, so to account for the positive part function the domain of integration is restricted to the area to the right of the dashed line (shaded area in the sketch).

The integral splits in 1D integrals in $dx_1$ and $dx_2$ as follows.

First it is the integral in $dx_1$, as shown by the strip in the sketch:

$ F(x_2) = int_{xi}^{infty} dx_1 rho_1(x_1) rho_2(x_2) (x_1 + x_2 - K) = F_1 + F_2 + F_3, $

where $xi = K - x_2$ for $x_2 < K$, and $xi=0$ for $x_2 > K$, and

$ F_1(x_2) = rho_2(x_2) int_{xi}^{infty} dx_1 rho_1(x_1) x_1 = rho_2(x_2) I_1(x_2), F_2(x_2) = x_2 rho_2(x_2) int_{xi}^{infty} dx_1 rho_1(x_1) = x_2 rho_2(x_2) I_2(x_2), F_3(x_2) = -K rho_2(x_2) int_{xi}^{infty} dx_1 rho_1(x_1) = -K rho_2(x_2) I_2(x_2), $

Once the functions $I_1$ and $I_2$ are known (i.e., evaluated to any desired accuracy) the second integral can be calculated as

$ int_0^{infty} dx_2 (F_1 + F_2 + F_3) $

Note that the integrals in $dx_1$, e.g., $I_{1}$ can be calculated as a running sum ($I_{1,0}=rho_{1,0}$, $I_{1,1}=I_{1,0}+rho_{1,1}$, $I_{1,k+1}=I_{1,k}+rho_{1,k+1}$, where $k$ is the grid index). That means that calculating it only once (using ${mathcal{O}}(N_1)$ operations) is enough to make $I_1(x_2)$ available for any $x_2$ by interpolation. For the 2D integral we'd need to use ${mathcal{O}}(N_2)$ operations on top of it, so the total count of mathematical operations is ${mathcal{O}}(N_1 + N_2) sim {mathcal{O}}(N)$, where $N$ is the characteristic grid size for one of the dimensions. If the integral was calculated directly as a 2D quadrature (i.e., a double sum) that would take a much larger number of mathematical operations, on the order of ${mathcal{O}}(N^2)$. It is not clear if the presented here approach could be extended beyond 2D; but if it could then the scaling of computational complexity would improve dramatically, from $mathcal{O}(N^M)$ to $mathcal{O}(Ntimes M)$, where $N$ is the characteristic grid size in one of the dimensions and $M$ is the number of dimensions.

Just another idea how you may proceed:

I restricted myself to 4 dimensions, however this should also work for 20 dimensions.

1. Use the Matlab symbolic Toolboox and define symbolic variables:

   syms x1 x2 x3 x4 ... k
   

2. Define your integration kernel:

   .

   $f(x_1, x_2, x_3, x_4)={mathrm{e}}^{-{x_{1}}^2},{mathrm{e}}^{-{x_{2}}^2},{mathrm{e}}^{-{x_{3}}^2},{mathrm{e}}^{-{x_{4}}^2},,left(k+x_{1}+x_{2}+x_{3}+x_{4}right)$

   .

   % Each summand separatly:
   
   f0(x1,x2,x3,x4)=exp(-x1^2)*exp(-x2^2)*exp(-x3^2)*exp(-x4^2)*k
   
   f1(x1,x2,x3,x4)=exp(-x1^2)*exp(-x2^2)*exp(-x3^2)*exp(-x4^2)*x1
   
   f2(x1,x2,x3,x4)=exp(-x1^2)*exp(-x2^2)*exp(-x3^2)*exp(-x4^2)*x2
   
   f3(x1,x2,x3,x4)=exp(-x1^2)*exp(-x2^2)*exp(-x3^2)*exp(-x4^2)*x3
   
   f4(x1,x2,x3,x4)=exp(-x1^2)*exp(-x2^2)*exp(-x3^2)*exp(-x4^2)*x4
   

3. Calculate the anti derivative:

   Looks straight forward, since $F$ is only a function of $boldsymbol{text{erf}}$ and $boldsymbol{text{exp}}$. You may define $F$ at once, however you won't see the pattern.

   % Integrate each summand separatly:
   
   F0(x1,x2,x3,x4)=int(int(int(int(f0,x1),x2),x3),x4)
   
                  =(k*pi^2*erf(x1)*erf(x2)*erf(x3)*erf(x4))/16
   
   F1(x1,x2,x3,x4)=int(int(int(int(f1,x1),x2),x3),x4) 
   
                  =-(pi^(3/2)*exp(-x1^2)*erf(x2)*erf(x3)*erf(x4))/16
   
   F2(x1,x2,x3,x4)=int(int(int(int(f2,x1),x2),x3),x4)
   
                  =-(pi^(3/2)*exp(-x2^2)*erf(x1)*erf(x3)*erf(x4))/16
   
   F3(x1,x2,x3,x4)=int(int(int(int(f3,x1),x2),x3),x4)
   
                  =-(pi^(3/2)*exp(-x3^2)*erf(x1)*erf(x2)*erf(x4))/16
   
   F4(x1,x2,x3,x4)=int(int(int(int(f4,x1),x2),x3),x4)
   
                  =-(pi^(3/2)*exp(-x4^2)*erf(x1)*erf(x2)*erf(x3))/16
   

4. Convert your symbolic expression to numerical expression (see MATLAB). Then you might use it directly in Matlab or in C++ or Fortran.

5. Use the divergence theorem

   $$begin{aligned}int_a^b f(t), dt & = F(b)-F(a) ... & = oint_{S} mathbf{F} cdot mathbf{n} , dS end{aligned}$$

6. Apply a Riemann sum or Trapezoidal rule for the multi-dimensional line integral with $aequiv 0,b equiv infty$.

   Advantages:

   • You only have to evaluate $F$ pointwise and add many $Delta F$
   • You won't have to save many numbers, only the anti derivative, $F$ and $Delta F$
   • You only apply a summation

Regards