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

Consider a potentially high-dimensional (say, $$N$$ up to 20) integral of the form $$\int_0^\infty \rho_1(x_1)\rho_2(x_2) \cdots \rho_N(x_N) \bigg(x_1+x_2+\cdots+x_N -K\bigg)^+ \, dx_1 \cdots dx_N.$$ where $$(z)^+ = \max(0,z)$$ is the positive part of the argument $$z\in \mathbb R$$, $$K \in \mathbb R_{>0}$$ and all $$\rho_m(x_m)>0$$.

I want to evaluate this integral numerically based on a representation of the product density functions $$\rho_m(x_m)$$ on a grid $$\{x_{m,k_m}\}$$, which basically leads to the term (integration weights set to one, for simplicity):

$$\sum_{k_1 \cdots k_N} \rho_{1,k_1} \cdots \rho_{N,k_N} \bigg(x_{1,k_1}+\cdots+x_{N,k_N} -K\bigg)^+$$

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.

I'd appreciate any comments on how to solve this integral numerically.

• What is K? What is the + superscript? – Maxim Umansky Jul 18 '20 at 16:01
• $(z)^+ = \max(z, 0)$ is the positive-part function and $K$ is a positive constant. – davidhigh Jul 18 '20 at 16:07
• Ok, so the integrand is rather simple, let's call it $f$. So if we consider a 2D case then it will be evaluation of a double sum of the form $\sum_{i=0,\infty} \sum_{j=0,\infty} \rho_{1,i} \rho_{2,j} f_{i,j}$. We know that to evaluate the integral we'll have to do the double sum. Then what is the issue? – Maxim Umansky Jul 18 '20 at 16:14
• That I can't do a 20-dimensional sum. Yet, if the "+" was missing, the integral would separate and I could evaluate it by 20 one-dimensional integrations. The question is basically, if there are clever trivks to restore the separability even with the "+", i.e. with the positive part function. – davidhigh Jul 18 '20 at 16:52
• Thanks for the like :-)! The rules around here are a bit strange. Actually, I'd rather discuss about equality in the world wide web. However, we can also discuss your question. Why do you think you can't apply a numerical quadrature in a 20 dimensional space? Even with the "+" part this should work in a tensorproduct like fashion. – ConvexHull Jul 18 '20 at 17:37

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:

$$\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 + 1\right)^2 \left[L_{n+1}\left(x_i\right)\right]^2}$$

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

$$\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?

• Thanks for your answer, and +1 so that you can comment in the future :-) Nevertheless, I try to avoid a direct numerical quadrature -- whatever the method -- as this won't work in the high-dimensional setting. – davidhigh Jul 18 '20 at 17:21

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}(N\times M)$$, where $$N$$ is the characteristic grid size in one of the dimensions and $$M$$ is the number of dimensions.

• Thanks for your answer. This approach amounts to direct integration, but with appropriate boundaries. In principle, what you sketch this is an application of the Fubini theorem, However, it scales as $\mathcal O(M^N)$ and is unfeasible in higher dimensions. – davidhigh Jul 18 '20 at 20:19
• Is that really true in this case? Do you agree that in 2D the number of function evaluations is $O(N_1 + N_2)$ rather than $O(N_1 \times N_2)$? – Maxim Umansky Jul 19 '20 at 4:14
• Thanks for your effort. In your edit, you neglected the positive part function $()^+$, and yes, in this case the effort boils down to $O(N_1+N_2)$. My question is now whether a similar thing is possible also with the "$+$". The direction I was thinking about (--as you asked about it in one former comment--) were integral transformations (Fourier, Laplace, and the like) where the coordinate-sum might factor, or other coordinate transformations (e.g. center-of-mass coordinates), or other clever tricks. – davidhigh Jul 19 '20 at 9:23
• I did not neglect the positive part function, the whole point of restricting the domain to the area to the right of the line $x_1 + x_2 - K = 0$ was to account for the positive part function. – Maxim Umansky Jul 19 '20 at 15:27
• Sorry for the delay. Your first formula amounts to an integration over $x_1$ for each $x_2$, thus yielding an $O(M^2)$ scheme. For larger dimensions, this function would become $F(x_2,…,x_N)$ and thus lead to an effort of $O(M^N)$ (where $M$=gridsize, $N$=dimension) – davidhigh Jul 20 '20 at 18:54

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


.

$$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 $$a\equiv 0,b \equiv \infty$$.

• 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$$