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.

  • $\begingroup$ What is K? What is the + superscript? $\endgroup$ Jul 18, 2020 at 16:01
  • $\begingroup$ $(z)^+ = \max(z, 0)$ is the positive-part function and $K$ is a positive constant. $\endgroup$
    – davidhigh
    Jul 18, 2020 at 16:07
  • $\begingroup$ 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? $\endgroup$ Jul 18, 2020 at 16:14
  • $\begingroup$ 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. $\endgroup$
    – davidhigh
    Jul 18, 2020 at 16:52
  • $\begingroup$ 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. $\endgroup$
    – ConvexHull
    Jul 18, 2020 at 17:37

3 Answers 3


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),$$


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

Hermite quadrature:

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


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


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?

  • 1
    $\begingroup$ 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. $\endgroup$
    – davidhigh
    Jul 18, 2020 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).

enter image description here

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.

  • $\begingroup$ 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. $\endgroup$
    – davidhigh
    Jul 18, 2020 at 20:19
  • $\begingroup$ 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)$? $\endgroup$ Jul 19, 2020 at 4:14
  • $\begingroup$ 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. $\endgroup$
    – davidhigh
    Jul 19, 2020 at 9:23
  • $\begingroup$ 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. $\endgroup$ Jul 19, 2020 at 15:27
  • $\begingroup$ 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) $\endgroup$
    – davidhigh
    Jul 20, 2020 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
  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:
  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:
  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$
    • You only apply a summation


  • $\begingroup$ Is this a real recipe or just random thoughts on the topic? If it's the first, could you elaborate on how to solve the integral? Thanks in advance. $\endgroup$
    – davidhigh
    Jul 20, 2020 at 18:45
  • $\begingroup$ It was a serious answer. However, you still have to apply a 19 dimensional quadrature of a hypercube surface. You only safe one dimension. $\endgroup$
    – ConvexHull
    Jul 21, 2020 at 5:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.