Local Versus Global Decisions in Bayesian Automatic Adaptive Quadrature

. The paper reports new signiﬁcant enhancement of the robustness and e ﬀ ective-ness of the Bayesian automatic adaptive quadrature over macroscopic integration ranges. The implementation of a classical m -panel rule (CC-32, Clenshaw-Curtis quadrature of algebraic degree of precision m = 32) is thought again. It involves new global and local decisions blocks which, on the one side, provide sharp diagnostics redirecting the advancement to the solution and, on the other side, take advantage of the progress in the available hardware to accelerate and to increase the accuracy of the computations. Where the decision power of CC-32 is exhausted, identiﬁcation and precise characterization of the features of the integrand proﬁle which prevent quick convergence are obtained by means of three-point Simpson rules spanned at triplets of successive CC-32 knots. This local complementary investigation tool provides scale insensitive diagnostics concerning the occurrence of integrand irregularities and prevents the activation of inappropriate decision blocks which would result in fake outputs.


Introduction
With a history starting at the time of Newton, the numerical solution of the one-dimensional Riemann integrals is the oldest, the most studied, and as yet only partially solved problem of the numerical analysis. Assuming that the integral of interest exists and has a finite value, I, the fundamental ingredient of the derivation of a numerical solution consists in the replacement of the original integrand by a simpler approximate expression for which the given integral can be solved analytically. The approximation of the integrand by an interpolatory polynomial is the most used approach to the derivation of a numerical solution. It results in an interpolatory quadrature sum, Q, which makes sense provided an associated error estimate, E > 0, is derived to assess its accuracy.
If the error estimate E exceeds the error tolerance specified at input, there are two possible approaches toward the improvement of the quality of the {Q, E > 0} pair: either by increasing the degree of the interpolatory polynomial over the input integration domain (see, [1,2] for more recent references) or by keeping a same quadrature sum and proceeding to the subdivision of the integration e-mail: adamg@jinr.ru,adamg@nipne.ro e-mail: adams@jinr.ru,adams@nipne.ro domain by successive subrange bisections (the m-panel rule approach either in the standard automatic adaptive quadrature (SAAQ), see, e.g., [3,4], or in the Bayesian automatic adaptive quadrature (BAAQ) over macroscopic ranges, see, e.g., [5] and references therein).
The present paper reports significant BAAQ robustness enhancement over macroscopic integration domains by supplementing the basic extended m-panel rule by a local quadrature rule.
For the m-panel rule of reference we choose the quadrature rule of high algebraic degree of precision, conveniently abridged as CC-32 (Clenshaw-Curtis quadrature of algebraic degree of precision m = 32) -for a motivation of this choice see [5]. While CC-32 accurately solves integrals involving smooth and moderately oscillatory integrands over moderately large integration domains, under the infringement of any these conditions, spurious quadrature sum Q outputs are obtained with no evidence for the reason of the bad output. The same is true for any other m-panel rule of high algebraic degree of precision, like the Gauss-Kronrod quadrature rules advocated in [3,4].
In the sequel, we propose to supplement the CC-32 rule with a sequence of local Simpson rules defined over triplets of successive knots of the CC-32 rule. The small (two-elementary-interval) extensions of the local Simpson rules enable sharp location and precise characterization of the irregular integrand profile features which ask for the redefinition of the path to the solution.
The paper is organized as follows. In section 2, basic BAAQ decision blocks are discussed. In section 3, generalized Simpson rules are derived over triplets of successive CC-32 knots. In section 4, selected empirical evidence is provided which illustrates the BAAQ robustness and SAAQ fragility. The paper ends with conclusions in section 5.

Basic decision blocks in the Bayesian automatic adaptive quadrature
We want to derive a BAAQ numerical solution of the (proper or improper) Riemann integral under the assumption that the real valued integrand function f (x) is continuous almost everywhere on the integration domain such that I exists and is finite. The weight function g(x) either absorbs an analytically integrable difficult factor in the integrand (e.g., endpoint singularity or oscillatory function) or else g( The automatic adaptive quadrature of (1) uses recurrence for the derivation of a global solution (Q N , E N > 0) satisfying the global heuristic termination criterion [3] where ε r and ε a denote the relative, respectively absolute, accuracy parameters requested at input. At the first step of the recurrence process, the m-panel rule yields a solution (Q 1 , E 1 > 0) of (1) over the whole integration domain [a, b]. The computation ends if this pair satisfies the termination criterion (2). Otherwise, subrange subdivision is performed with the aim at improving the accuracy of the numerical solution. This generates a partition of [a, b], and over each subrange of this partition, the m-panel rule is used to compute a local numerical solution (q k , e k > 0). The global solution of (1) is obtained by summing up the local solutions, In the expression of Q N , q 0 and q N+1 denote input values defined as follows: q 0 = 0 if a > −∞ and q N+1 = 0 if b < +∞. Otherwise, in the partition (3), x 0 > a = −∞ and x N < b = +∞ denote finite cut-offs of (1) provided as inputs to the numerical problem together with q 0 0 and q N+1 0 denoting integrand dependent asymptotic contributions to Q N over (−∞, x 0 ] and [x N , +∞), respectively.
The SAAQ is based on two fundamental building blocks. First, an m-panel rule of high algebraic degree of precision is used for the approximation of the integral (1) irrespective of the extension and location of the integration domain [a, b] over the real axis. Second, the refinement of the partition (3) is achieved by subrange bisection such that the generated subranges define a complete binary tree the evolution of which is controlled by an associated priority queue.
It was noticed long time ago [6] that this m-panel automatic adaptive quadrature (AAQ) approach to the solution may fail badly. To circumvent the existing difficulties, the SAAQ solution proposed in QUADPACK [3] was to implement separate AAQ codes for distinct classes of Riemann integrals (1) leaving to the user the responsibility to choose the convenient code from the proposed menu. The QUADPACK SAAQ was implemented in the major numerical software packages like NAG [7]. A step forward was made in the open source Maxima symbolic-numerical package [8]. Here the user inputs the problem through a convenient graphical user interface (GUI) which selects automatically the most convenient code from the existing SLATEC implementation of QUADPACK [9].
Within the BAAQ, as documented elsewhere [10,11], under floating point computations, the use of an m-panel rule of high algebraic degree of precision over a subrange [α (3)  To evaluate its magnitude, entailing either EOC (under catastrophic cancellation by subtraction) or adjustments of the input accuracy parameters ε a and ε r to the (lower) ceiling of the available accuracy, the driver procedure invokes an analysis decision block using a composite Simpson rule with centroid spanned local Simpson rules implemented at the Fejer subset [5] of the CC-32 rule (Sec. 3).
I.4. If no EOC, computation of the Chebyshev expansion coefficients of the CC-32 rule is done by a fast Chebyshev transform (FCT) [12], together with their conditioning diagnostic.
I. 5. Under the return of a well-conditioning diagnostic of the FCT output, activation of the CC-32 rule over [a, b] and check of the termination criterion (2) for (Q 1 , E 1 > 0) are done.
II. Decision blocks on subrange subdivisions are activated whenever the Bayesian diagnostics ask for such operations or the termination criterion (2) is not fulfilled. This is not the subdivision by bisection of the SAAQ. The previously acquired knowledge on the IP generated over the parent subrange is used to predict the most appropriate partition into as many descendent subranges as necessary.
The inheritance property: the ends of each descendent subrange are set at conveniently defined knots of the parent subrange. The reduced length of a descendent is known to processor accuracy. It is inherited together with all the relevant already computed quantities at the subrange ends.
II.1. If a piecewise linear IP is resolved (the signature of which is the occurrence of pieces of three or more consecutively vanishing local IP curvatures), then immediate subrange subdivision is performed. This results in alternating linear and nonlinear descendent subranges. The conditioning of the nonlinear descendents is subject to further study. [13][14][15]. As a result of the arcsine distribution of the quadrature knots of the CC-32 rule [2], the highest knot densities occur inside the CPNSEs which provide thus the most reliable information on the IP variation. Each of the two independently done CPNSE analyses has two main targets.

II.2. Inspections of the close proximity neighbourhoods of the integration domain ends (CPNSEs)
(i) Assessment of the integrand conditioning at each of the two endpoints of the integration domain. The ratios of the local slopes computed over the nearest neighbouring and the next nearest neighbouring elementary intervals with respect to a same endpoint yield two possible diagnostics: Rregular behaviour (ratio close to unity); IRirregular behaviour (ratio much larger or much smaller than one). The latter diagnostic needs refinement to discriminate between the following three cases: Bboundary layer, i.e., singularity at the endpoint of interest (asks for appropriate activation of a convergence acceleration procedure); Ffake boundary layer, regular behaviour occurs over smaller integration intervals (selects the subrange subdivision as the only robust path to the derivation of the solution); Jendpoint jump of zero measure. This refinement asks for a few (two or three) additional integrand evaluations, inside the integration domain, at relative distances equating the machine epsilon with respect to the addition and comparison of the integrand variations over such equally spaced tiny distances [16].
(ii) Assessment of the integrand conditioning over each inner CPNSE region from the behaviour pattern of local slopes and local curvatures returns one of the following three diagnostics: Rregular nonlinear IP behaviour; Oviolent integrand oscillations the characteristic signature of which is the presence of two IP extrema at least inside a CPNSE; SII -characteristic signature of an isolated irregularity (inner singularity, two-sided asymptote, turning point, or finite jump).
The combinations RO, OR, or OO ask for the use of a subrange subdivision procedure of the current parent range which shares the number of already resolved oscillations over descendents in (3).
A combined RR diagnostic may point to an inner irregular integrand behaviour. It is to be stressed the tentative character of the SII Bayesian diagnostics. If the iteration of the CC-32 rule over such a subrange results in further enhancement of the amplitude of the irregular variation pattern, the decision is taken to locate, as fast as possible, the irregular feature to processor accuracy and to use an adequate convergence acceleration procedure afterwards. If attenuation or disappearance of the irregular variation pattern occurs after the CC-32 rule iteration, the pattern for solving regular integrands is followed over the descendents of the corresponding subrange. Sampling of integration interval "Ivals1" "S1-s10" "C1-s100"

Three-point Simpson rules over CC-32 quadrature knots
In this section, three-point Simpson rules of local third order accuracy are derived over triplets of successive (unequally spaced) knots of the CC-32 rule.
Since the Simpson rules have local spanning support, their addition to the CC-32 rule enables the identification and precise characterization of those local CC-32 IP features which prevent reliable and accurate CC-32 solution. The Facts 3.1 and 3.2 stated below direct the formulation of the three-point Simpson rules such as to secure the maximum available floating point accuracy. With the original convention used by Clenshaw and Curtis [17], we have the set of standard reduced CC quadrature knots (SRKs), y 4n j = cos(π j/4n), j = 0, 1, . . . , 4n. The floating point compu-tations of the local slopes and curvatures entail precision losses in the denominators of the first and second order divided differences, which get maximized inside the CPNSEs of the ends α and β. The alternative use of the MRKs [5] avoids this annoying precision loss in the computation of the divided differences. The MRKs are defined as the absolute distances from the SRKs to the nearest endpoints of the reduced range [−1, 1] (inside (0, 1], η 2n j = 1 − y 4n j , j = 0, 1, . . . , 2n − 1; inside [−1, 0), η 2n j = 1 + y 4n 4n− j , j = 0, 1, . . . , 2n − 1; at zero: η 2n 2n = 1). The MRKs over each half range are given by 0 = η 2n 0 < η 2n 1 < · · · η 2n j < · · · < η 2n 2n = 1.
Over the triplet centred at η 2n 2n = 1, the second term in Eq. (9) vanishes and we recover the standard three point Simpson rule over equally spaced knots.
While algebraically equivalent, the lCS, Eq. (9), and lTS, Eq. (10), rules show different numerical conditioning properties. In view of (7), all the coefficients of the centroid spanned three-point Simpson rule (9) are positive. Therefore, this rule may be confidently used for the assessment of the magnitude of the cancellation by subtraction discussed in Sec. 2. On the other side, in view of the additivity property of the trapezoidal rule, the trapezoidal spanned three-point Simpson rule (10) can be used for the derivation of quadrature error estimates.
II. Local quadrature error estimates are obtained by elementary operations. For instance, the three-point Simpson rule over the interval [η 2n j−1 , η 2n j+1 ] associates the local error estimate

BAAQ is robust, SAAQ is fragile
Perhaps the most remarkable feature of practical interest of the present BAAQ implementation is its robustness to unsuitable formulation of the integral (1). Two relevant case studies are discussed. Case Study 1. Bringing an isolated inner integrand singularity to subrange integration ends. [3] (p. 110), shows an inner integrand singularity at x s = √ 3 − 1 0.7320. The splitting of the integral over [0, x s ] and [x s , 1] yields, by means of the DQAGS subroutine of QUADPACK (and its corresponding NAG implementation [7]), reliable output due to the appropriate activation of the epsilon acceleration algorithm [3], Sec. 2.2.4. However, if the singularity was not given as an input to the code, the bisection based subrange subdivision eventually results in a low accuracy output together with a pointer diagnostic showing that an unreliable solution was obtained.
The wxMaxima automatically selects the quad_qags function with results similar to QUADPACK. The workaround proposed in QUADPACK for this case ( [3], pp. 110-111) is to activate the subroutine DQAGE which yields, after 435 integrand evaluations, the message that a local difficulty occurs in the interval (0.73199, 0.73206). This returns the problem to the user for separate analysis.
Within the present BAAQ implementation, a three-step fast singularity automatic location is used: (i) Bayesian inference (CC-32 solution over [0, 1] yields irregular behavior of Chebyshev expansion coefficients. Activation of three-point Simpson rules finds a singularity signature within the interval (0.71, 0.75) (Fig. 1 -left)); (ii) inference confirmation (iteration of the procedure over this interval confirms the singularity signature at some abscissa x s within a narrower subinterval of approximate width 0.0003 ( Fig. 1 -right)); (iii) resolving the value x s to machine accuracy, using an extremum search procedure. Afterwards, the automatic splitting of the integral as above yields accurate solution.
Case Study 2. Identification of a fake boundary layer prevents inappropriate activation of the convergence acceleration procedure.
The integral (1) with a = 0, b = ∞, f (x) = (x 2 + 1) −1 (g(x) = 1) [18] has the value π/2. A rapid look at the integral shows that a finite cut-off T can be defined such that, at asymptotically large x, T < x < ∞, the unity can be neglected with respect to x 2 . The integral over [T, ∞) yields the asymptotic contribution q N+1 = 1/T to the expression (4) of Q N . If machine accuracy of the numerical solution is sought, then the processor accuracy defined threshold T = 2 32 is to be used.
The SAAQ solution to the finite integral over [0, T ], T 1, was shown several decades ago [18] to yield a spurious negative value with the wrong diagnostic of reliable output due to the inappropriate decision to activate the convergence acceleration procedure.
The wxMaxima solution of this integral over [0, +∞) activates the function quad_qagi which yields correct output. However, the attempt at solving the integral over [0, T ] by means of wxMaxima activates the function quad_qags and fails.
The present local-Simpson-rule based analysis identifies a fake boundary layer at the origin, hence the ban to activate the convergence acceleration procedure. Then a numerical solution based on mere subrange subdivision results in the correct output of π/2 to within machine accuracy.
Due to the dependence of the SAAQ output on the parameters assigned to the formulation of the Riemann integral, we conclude that the SAAQ approach to the numerical solution of the Riemann integral (1) is fragile. In contradistinction to it, the present BAAQ implementation yields reliable accurate outputs for the integral (1). The BAAQ algorithm is robust.

Conclusions
The discussed Bayesian automatic adaptive quadrature (BAAQ) over integration ranges of macroscopic extension bears little resemblance to the standard automatic adaptive quadrature (SAAQ).
Within the BAAQ, both the understanding of the importance of the different features of the integrand profile (IP) and the evolution of the hardware at hand have resulted in new organization of the IP computation and decisions on the iterative advancement to the solution. Moreover, the combination of the m-panel rule with a conditionally activated rule encompassing local three-point Simpson rules defined over triplets of adjacent knots of the reference m-panel rule enables identification and precise characterization of the IP features which prevent the derivation of the reliable output.
In contradistinction to the SAAQ, the derivation of an m-panel rule output over the whole integration domain is subject to a complex Bayesian decision process which gradually unveils the essential Riemann integral features during the computations and redirects the computing chain accordingly.
If subrange subdivision is needed, the use of the parameters of the local Simpson rules enables fine tuning of the subrange subdivision to the particularities of the variation of the integrand function such that, both for signatures pointing to irregular IP variations or to regular ones, robust decisions are taken which result in reliable outputs under considerable acceleration of the computing speed.