Introduction
In this document are given outlines and examples of several related implementations of Lebesgue integration, [1], within the framework of NIntegrate
, [7]. The focus is on the implementations of Lebesgue integration algorithms that have multiple options and can be easily extended (in order to do further research, optimization, etc.) In terms of NIntegrate
‘s framework terminology it is shown how to implement an integration strategy or integration rule based on the theory of the Lebesgue integral. The full implementation of those strategy and rules — LebesgueIntegration
, LebesgueIntegrationRule
, and GridLebesgueIntegrationRule
— are given in the Mathematica package [5].
The advantage of using NIntegrate
‘s framework is that a host of supporting algorithms can be employed for preprocessing, execution, experimentation, and testing (correctness, comparison, and profiling.)
Here is a brief description of the integration strategy LebesgueIntegration
in [5]:
 prepare a function that calculates measure estimates based on random points or low discrepancy sequences of points in the integration domain;

use
NIntegrate
for the computation of one dimensional integrals for that measure estimate function over the range of the integrand function values.
The strategy is adaptive because of the second step — NIntegrate
uses adaptive integration algorithms.
Instead of using an integration strategy we can "tuck in" the whole Lebesgue integration process into an integration rule, and then use that integration rule with the adaptive integration algorithms NIntegrate
already has. This is done with the implementations of the integration rules LebesgueIntegrationRule
and GridLebesgueIntegrationRule
.
Brief theory
Lebesgue integration extends the definition of integral to a much larger class of functions than the class of Riemann integrable functions. The Riemann integral is constructed by partitioning the integrand’s domain (on the axis). The Lebesgue integral is constructed by partitioning the integrand’s codomain (on the axis). For each value in the codomain, find the measure of the corresponding set of points in the domain. Roughly speaking, the Lebesgue integral is then the sum of all the products ; see [1]. For our implementation purposes is defined differently, and in the rest of this section we follow [3].
Consider the nonnegative boundable measurable function :
We denote by the measure for the points in for which , i.e.
The Lebesgue integral of over can be be defined as:
Further, we can write the last formula as
The restriction can be handled by defining the following functions and :
and using the formula
Since finding analytical expressions of is hard we are going to look into ways of approximating .
For more details see [1,2,3,4].
(Note, that the theoretical outline the algorithms considered can be seen as algorithms that reduce multidimensional integration to one dimensional integration.)
Algorithm walk through
We can see that because of Equation (4) we mostly have to focus on estimating the measure function . This section provides a walk through with visual examples of a couple of stochastic ways to do that.
Consider the integral
In order to estimate in for we are going to generate in a set of low discrepancy sequence of points, [2]. Here this is done with points of the so called "Sobol" sequence:
n = 100;
SeedRandom[0,Method > {"MKL",Method > {"Sobol", "Dimension" > 2}}];
points = RandomReal[{0, 1}, {n, 2}];
ListPlot[points, AspectRatio > Automatic,
PlotTheme > "Detailed",
FrameLabel > {"\!\(\*SubscriptBox[\(x\), \(1\)]\)","\!\(\*SubscriptBox[\(x\), \(2\)]\)"}]
To each point let us assign a corresponding "volume" that can be used to approximate with Equation (2). We can of course easily assign such volumes to be , but as it can be seen on the plot this would be a good approximation for a larger number of points. Here is an example of a different volume assignment using a Voronoi diagram, [10]:
vmesh = VoronoiMesh[points, {{0, 1}, {0, 1}}, Frame > True];
Show[{vmesh, Graphics[{Red, Point[points]}]},
FrameLabel > {"\!\(\*SubscriptBox[\(x\), \(1\)]\)", "\!\(\*SubscriptBox[\(x\), \(2\)]\)"}]
(The Voronoi diagram for finds for each point the set of domain points closest to than any other point of .)
Here is a breakdown of the Voronoi diagram volumes corresponding to the generated points (compare with ) :
volumes =PropertyValue[{vmesh, Dimensions[points][[2]]}, MeshCellMeasure];
Histogram[volumes, PlotRange > All, PlotTheme > "Detailed",
FrameLabel > {"volume", "count"}]
Let us define a function that computes according to Equation (1) with the generated points and assigned volumes:
EstimateMeasure[fval_?NumericQ, pointVals_, pointVolumes_] :=
Block[{pinds},
pinds = Clip[Sign[pointVals  fval], {0, 1}, {0, 1}];
pointVolumes.pinds
];
Here is an example call of that function using the Voronoi diagram volumes:
EstimateMeasure[1.6, Sqrt[2 + Total[#]] & /@ points, volumes]
(* 0.845833 *)
And here is another call using uniform volumes:
EstimateMeasure[1.6, Sqrt[2 + Total[#]] & /@ points,
Table[1/Length[points], {Length[points]}]] // N
(* 0.85 *)
The results can be verified using ImplicitRegion
:
RegionMeasure[ImplicitRegion[ Sqrt[2 + x1 + x2] >= 1.6, {{x1, 0, 1}, {x2, 0, 1}}]]
(* 0.8432 *)
Or using Integrate
:
Integrate[Piecewise[{{1, Sqrt[2 + x1 + x2] >= 1.6}}, 0], {x1, 0, 1}, {x2, 0, 1}]
(* 0.8432 *)
At this point we are ready to compute the integral estimate using Formula (4) :
fvals = Sqrt[2 + Total[#]] & /@ points;
Min[fvals]*EstimateMeasure[0, fvals, volumes] +
NIntegrate[EstimateMeasure[y, fvals, volumes], {y, Min[fvals], Max[fvals]}]
(* 1.72724 *)
To simplify the code we use the symbol to hold the values . Note that instead of the true min and max values of we use estimates of them with .
Here is the verification of the result:
Integrate[Sqrt[2 + x1 + x2], {x1, 0, 1}, {x2, 0, 1}]
% // N
(* 8/15 (16 + 2 Sqrt[2]  9 Sqrt[3]) *)
(* 1.72798 *)
In order to implement the outlined algorithm so it will be more universal we have to consider volumes rescaling, function positivity, and Voronoi diagram implementation(s). For details how these considerations are resolved see the code of the strategy LebesgueIntegration in [5].
Further algorithm elaborations
Measure estimate with regular grid cells
The article 3 and book 4 suggest the measure estimation to be done through membership of regular grid of cells. For example, the points generated in the previous section can be grouped by a grid:
The following steps describe in detail an algorithm based on the proposed in [3,4] measure estimation method. The algorithm is implemented in [5] for the symbol, GridLebesgueIntegrationRule
.
1. Generate points filling the where is the dimension of .
2. Partition with a regular grid according to specifications.
 2.1. Assume the cells are indexed with the integers , .
 2.2. Assume that all cells have the same volume below.
3. For each point find to which cell of the regular grid it belongs to.
4. For each cell have a list of indices corresponding to the points that belong to it.
5. For a given subregion of integration rescale to the points to lie within ; denote those points with .
 5.1. Compute the rescaling factor for the integration rule; denote with .
6. For a given integrand function evaluate over all points .
7. For each cell find the min and max values of .
 7.1. Denote with and correspondingly.
8. For a given value , where is some integer enumerating the 1D integration rule sampling points, find the coefficients , using the following formula:
9. Find the measure estimate of with
Axis splitting in Lebesgue integration rules
The implementations of Lebesgue integration rules are required to provide a splitting axis for use of the adaptive algorithms. Of course we can assign a random splitting axis, but that might lead often to slower computations. One way to provide splitting axis is to choose the axis that minimizes the sum of the variances of subdivided regions estimated by samples of the rule points. This is the same approach taken in NIntegrate`s rule "MonteCarloRule"; for theoretical details see the chapter "7.8 Adaptive and Recursive Monte Carlo Methods" of [11].
In [5] this splitting axis choice is implemented based on integrand function values. Another approach, more in the spirit of the Lebesgue integration, is to select the splitting axis based on variances of the measure function estimates.
Consider the function:
DensityPlot[Exp[3 (x  1)^2  4 (y  1)^2], {x, 0, 3}, {y, 0, 3},
PlotRange > All, ColorFunction > "TemperatureMap"]
Here is an example of sampling points traces using "minimum variance" axis splitting in LebesgueIntegrationRule
:
res = Reap@
NIntegrate[Exp[3 (x  1)^2  4 (y  1)^2], {x, 0, 3}, {y, 0, 3},
Method > {"GlobalAdaptive",
Method > {LebesgueIntegrationRule, "Points" > 600,
"PointGenerator" > "Sobol",
"AxisSelector" > "MinVariance"},
"SingularityHandler" > None},
EvaluationMonitor :> Sow[{x, y}],
PrecisionGoal > 2.5, MaxRecursion > 3];
res[[1]]
ListPlot[res[[2, 1]], AspectRatio > Automatic]
(* 0.890916 *)
And here are the sampling points with random selection of a splitting axis:
res = Reap@
NIntegrate[Exp[3 (x  1)^2  4 (y  1)^2], {x, 0, 3}, {y, 0, 3},
Method > {"GlobalAdaptive",
Method > {LebesgueIntegrationRule, "Points" > 600,
"PointGenerator" > "Sobol",
"AxisSelector" > Random},
"SingularityHandler" > None},
EvaluationMonitor :> Sow[{x, y}],
PrecisionGoal > 2.5, MaxRecursion > 3];
res[[1]]
ListPlot[res[[2, 1]], AspectRatio > Automatic]
(* 0.892499 *)
Here is a more precise estimate of that integral:
NIntegrate[Exp[3 (x  1)^2  4 (y  1)^2], {x, 0, 3}, {y, 0, 3}]
(* 0.898306 *)
Implementation design within NIntegrate’s framework
Basic usages
The strategy and rule implementations in [5] can be used in the following ways.
NIntegrate[Sqrt[x], {x, 0, 2}, Method > LebesgueIntegration]
(* 1.88589 *)
NIntegrate[Sqrt[x], {x, 0, 2},
Method > {LebesgueIntegration, Method > "LocalAdaptive",
"Points" > 2000, "PointGenerator" > "Sobol"}, PrecisionGoal > 3]
(* 1.88597 *)
NIntegrate[Sqrt[x], {x, 0, 2},
Method > {LebesgueIntegrationRule, "Points" > 2000,
"PointGenerator" > "Sobol",
"PointwiseMeasure" > "VoronoiMesh"}, PrecisionGoal > 3]
(* 1.88597 *)
NIntegrate[Sqrt[x + y + x], {x, 0, 2}, {y, 0, 3}, {z, 0, 4},
Method > {GridLebesgueIntegrationRule,
Method > "GaussKronrodRule", "Points" > 2000,
"GridSizes" > 5, "PointGenerator" > "Sobol"}, PrecisionGoal > 3]
(* 43.6364 *)
Options
Here are the options for the implemented strategy and rules in [5]:
Grid[Transpose[{#, ColumnForm[Options[#]]} & /@
{LebesgueIntegration,LebesgueIntegrationRule,
GridLebesgueIntegrationRule}],
Alignment > Left, Dividers > All]
Using variable ranges
Integration with variable ranges works "out of the box."
NIntegrate[Sqrt[x + y], {x, 0, 2}, {y, 0, x}]
(* 2.75817 *)
NIntegrate[Sqrt[x + y], {x, 0, 2}, {y, 0, x},
Method > LebesgueIntegration, PrecisionGoal > 2]
(* 2.75709 *)
NIntegrate[Sqrt[x + y], {x, 0, 2}, {y, 0, x},
Method > LebesgueIntegrationRule, PrecisionGoal > 2]
(* 2.75663 *)
Infinite ranges
With infinite ranges the implemented Lebesgue integration algorithms produce estimates far from the real values. Here is an example:
NIntegrate[1/x^2, {x, 1, Infinity}]
(* 1. *)
res = Reap@
NIntegrate[1/x^2, {x, 1, \[Infinity]},
Method > {LebesgueIntegration, "Points" > 400,
"PointGenerator" > "Sobol",
"LebesgueIntegralVariableSymbol" > fval},
EvaluationMonitor :> {Sow[fval]},
PrecisionGoal > 3, MaxRecursion > 20,
WorkingPrecision > 20];
res = DeleteCases[res, fval, \[Infinity]];
res[[1]]
ListPlot[res[[2, 1]]]
(* 399.91652949396108627 *)
In order to achieve better results dedicated implementation changes (or separate development) have to be done.
Evaluation monitoring
With the option "EvaluationMonitor" we can see the sampling points for the strategy and the rules.
This is straightforward for the rules:
res = Reap@
NIntegrate[Exp[3 (x  1)^2], {x, 0, 3},
Method > LebesgueIntegrationRule,
EvaluationMonitor :> Sow[x],
PrecisionGoal > 3];
ListPlot[res[[2, 1]]]
The strategy LebesgueIntegration
uses an internal variable for the calculation of the Lebesgue integral. In "EvaluationMonitor" either that variable has to be used, or a symbol name has to be passed though the option "LebesgueIntegralVariableSymbol". Here is an example:
res = Reap@
NIntegrate[Sqrt[x + y + z], {x, 1, 2}, {y, 0, 1}, {z, 1, 12},
Method > {LebesgueIntegration, "Points" > 600,
"PointGenerator" > "Sobol",
"LebesgueIntegralVariableSymbol" > fval},
EvaluationMonitor :> {Sow[fval]},
PrecisionGoal > 3];
res = DeleteCases[res, fval, \[Infinity]];
ListPlot[res[[2, 1]]]
Profiling
We can use NIntegrate
‘s utility functions for visualization and profiling in order to do comparison of the implemented algorithms with related ones (like "AdaptiveMonteCarlo") which NIntegrate
has (or are pluggedin).
Needs["Integration`NIntegrateUtilities`"]
Common function and domain:
fexpr = 1/(375/100  Cos[x]  Cos[y]);
ranges = {{x, 0, \[Pi]}, {y, 0, \[Pi]}};
Common parameters:
pgen = "Random";
npoints = 1000;
"AdaptiveMonteCarlo" call:
NIntegrateProfile[
NIntegrate[fexpr, Evaluate[Sequence @@ ranges],
Method > {"AdaptiveMonteCarlo",
Method > {"MonteCarloRule", "PointGenerator" > pgen,
"Points" > npoints,
"AxisSelector" > "MinVariance"}},
PrecisionGoal > 3]
]
(* {"IntegralEstimate" > InputForm[2.8527356472097902`], "Evaluations" > 17000, "Timing" > 0.0192079} *)
LebesgueIntegrationRule
call:
NIntegrateProfile[
NIntegrate[fexpr, Evaluate[Sequence @@ ranges],
Method > {"GlobalAdaptive",
Method > {LebesgueIntegrationRule,
"PointGenerator" > pgen, "Points" > npoints,
"AxisSelector" > "MinVariance",
Method > "ClenshawCurtisRule"},
"SingularityHandler" > None}, PrecisionGoal > 3]
]
(* {"IntegralEstimate" > InputForm[2.836659588960318], "Evaluations" > 13000, "Timing" > 0.384246} *)
Design diagrams
The flow charts below show that the plugin designs have common elements. In order to make the computations more effective the rule initialization prepares the data that is used in all rule invocations. For the strategy the initialization can be much lighter since the strategy algorithm is executed only once.
In the flow charts the double line boxes designate subroutines. We can see that so called Hollywood principle "don’t call us, we’ll call you" in Objectoriented programming is manifested.
The following flow chart shows the steps of NIntegrate
‘s execution when the integration strategy LebesgueIntegration is used.
The following flow chart shows the steps of NIntegrate
‘s execution when the integration rule LebesgueIntegrationRule
is used.
Algorithms versions and options
There are multiple architectural, coding, and interface decisions to make while doing implementations like the ones in [5] and described in this document. The following mind map provides an overview of alternatives and interactions between components and parameters.
Alternative implementations
In many ways using the Lebesgue integration rule with the adaptive algorithms is similar to using NIntegrate
‘s "AdaptiveMonteCarlo" and its rule "MonteCarloRule". Although it is natural to consider pluggingin the Lebesgue integration rules into "AdaptiveMonteCarlo" at this point NIntegrate
‘s framework does not allow "AdaptiveMonteCarlo" the use of a rule different than "MonteCarloRule".
We can consider using Monte Carlo algorithms for estimating the measures corresponding to a vector of values (that come from a 1D integration rule). This can be easily done, but it is not that effective because of the way NIntegrate
handles vector integrands and because of stopping criteria issues when the measures are close to .
Future plans
One of the most interesting extensions of the described Lebesgue integration algorithms and implementation designs is their extension with more advanced features of Mathematica for geometric computation. (Like the functions VoronoiMesh
, RegionMeasure
, and ImplicitRegion
used above.)
Another interesting direction is the derivation and use of symbolic expressions for the measure functions. (Hybrid symbolic and numerical algorithms can be obtained as NIntegrate
‘s handling of piecewise functions or the strategy combining symbolic and numerical integration described in [9].)
References
[1] Wikipedia entry, Lebesgue integration, URL: https://en.wikipedia.org/wiki/Lebesgue_integration .
[2] Wikipedia entry, Lowdiscrepancy sequence, URL: https://en.wikipedia.org/wiki/Lowdiscrepancy_sequence .
[3] B. L. Burrows, A new approach to numerical integration, 1. Inst. Math. Applics., 26(1980), 151173.
[4] T. He, Dimensionality Reducing Expansion of Multivariate Integration, 2001, Birkhauser Boston. ISBN13:9781461274148 .
[5] A. Antonov, Adaptive Numerical Lebesgue Integration Mathematica Package, 2016, MathematicaForPrediction project at GitHub.
[6] A. Antonov, Lebesgue integration, Wolfram Demonstrations Project, 2007. URL: http://demonstrations.wolfram.com/LebesgueIntegration .
[7] "Advanced Numerical Integration in the Wolfram Language", Wolfram Research Inc. URL: https://reference.wolfram.com/language/tutorial/NIntegrateOverview.html .
[8] A. Antonov, "How to implement custom integration rules for use by NIntegrate?", (2016) Mathematica StackExchange answer, URL: http://mathematica.stackexchange.com/a/118326/34008 .
[9] A. Antonov, "How to implement custom NIntegrate integration strategies?", (2016) Mathematica StackExchange answer, URL: http://mathematica.stackexchange.com/a/118922/34008 .
[10] Wikipedia entry, Voronoi diagram, URL: https://en.wikipedia.org/wiki/Voronoi_diagram .
[11] Press, W.H. et al., Numerical Recipes in C (2nd Ed.): The Art of Scientific Computing, 1992, Cambridge University Press, New York, NY, USA.