Disentangling Complexity in Bayesian Automatic Adaptive Quadrature

. The paper describes a Bayesian automatic adaptive quadrature (BAAQ) so-lution for numerical integration which is simultaneously robust, reliable, and e ﬃ cient. Detailed discussion is provided of three main factors which contribute to the enhancement of these features: (1) reﬁnement of the m -panel automatic adaptive scheme through the use of integration-domain-length-scale-adapted quadrature sums; (2) fast early problem complexity assessment – enables the non-transitive choice among three execution paths: (i) immediate termination (exceptional cases); (ii) pessimistic – involves time and resource consuming Bayesian inference resulting in radical reformulation of the problem to be solved; (iii) optimistic – asks exclusively for subrange subdivision by bisection; (3) use of the weaker accuracy target from the two possible ones (the input accuracy speciﬁcations and the intrinsic integrand properties respectively) – results in maximum possible solution accuracy under minimum possible computing time.


Introduction
The present paper reports results along the lines of our recent investigations [1,2] concerning the Bayesian automatic adaptive quadrature (BAAQ) solution [3] of definite Riemann integrals.
The essential requirements to the BAAQ algorithms are robustness, reliability, and efficiency under floating point computations. Their implementation asks for three main BAAQ developments.
(a) The refinement of the m-panel scheme of the standard adaptive quadrature (described, e.g., in [4][5][6][7]). This development was suggested by the fact that, under floating point computations, the algebraic degree of precision of an interpolatory quadrature sum ceases to be a characteristic invariant feature of it. The floating point degree of precision is the adequate quantity characterizing a quadrature sum [8]. A convenient way of taking into account this floating point feature was [1] to modify the m-panel scheme by the definition of three classes of integration domain length scales: macroscopic, mesoscopic, and microscopic and to use characteristic quadrature sums over each scale.
(b) The fast early problem complexity assessment is a new undertaking the need of which follows from the empirical fact that a single algorithm cannot efficiently cope with the myriad of existing e-mail: adamg@jinr.ru,adamg@theory.nipne.ro e-mail: adams@jinr.ru,adams@theory.nipne.ro Riemann integrals. However, instead of leaving with the user the expert decision on the choice among the large number of existing computing codes for numerical solution of the integrals (see, [4][5][6][7]), an attempt was made to investigate the possibility to define decision paths based on Bayesian inferences which are tailoring automatically the best way to follow to the derivation of accurate output.
(c) Finally, new results are reported concerning the code optimization in connection with the user accuracy requirements. A useful by-product of the problem complexity assessment is the possibility to define an integrand adapted upper accuracy target. Whenever this is less stringent than the target following from the user defined accuracy requirements, the latter is superseded by the former one. This decision has two desirable effects: it provides at output the best possible accuracy for the problem at hand, while sparing the unnecessary computations done beyond the maximum available accuracy.
With this general plan in mind, the discussion below is organized in four sections. Section 2 provides a summary of the BAAQ approach enabling independent lecture of this paper. Section 3 adds progress to the three class m-panel scheme proposed in [1]. Section 4 describes in detail the complexity assessment over macroscopic ranges. It also summarizes the integrand induced accuracy limitations. The paper ends with concluding remarks in section 5.

Summary of the Bayesian automatic adaptive quadrature
A BAAQ numerical solution of the (proper or improper) Riemann integral is sought under the assumption that the real valued integrand function f (x) is continuous almost everywhere on [a, b] such that (1) exists and is finite. Here 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(x) ≡ 1, ∀x ∈ [a, b]. Following [1], the subsequent BAAQ developments implement a solution of (1) within a class dependent m-panel rule approach. We remind that the m-panel rule idea was proposed within the standard automatic adaptive quadrature (AAQ) (see, e.g., [4][5][6]) as a means to adjust the discretization of the original integration domain [a, b] As pointed out in [1], the quadrature sum q is to be computed by different quadrature sums over three classes of integration domain widths: macroscopic, mesoscopic, and microscopic, respectively.
An AAQ algorithm works iteratively. The quadrature rule at hand provides a global output initialization {Q N , E N |N = 1} over [a, b]. If the obtained approximation Q 1 of (1) is not accurate enough, then global output refinement is obtained by gradual subdivision of  [4], where, ε a and ε r denote the absolute accuracy, respectively relative accuracy requested at input.
The termination of the iterative process may end successfully, or it may not. The reason for failure is the evaluation of local error estimates by use of probabilistic arguments which might result in one or more persistently spurious e i values.
The BAAQ advancement to the solution rests on Bayesian inference based on four pillars: the theory of the Riemann integral, the theory of the numerical integration (quadrature), features of the floating point computation, the accumulated empirical evidence. Essentially, the probabilistic character of the AAQ approach to the derivation of local error estimates e > 0 is preserved. However, each step of the gradual advancement to the solution is scrutinized based on a set of hierarchically ordered criteria which enable decision taking in terms of the established diagnostics [3].

Refinement of the m-panel AAQ scheme
The inadequacy of the algebraic degree of precision for the characterization of the outputs of a quadrature sum in floating point computations [8] is illustrated in the figure 1. Extensive numerical tests have confirmed the usefulness of the m-panel AAQ scheme refinement proposed in [1], concerning the definition of three length scales of the integration domain widths: macroscopic, mesoscopic, and microscopic, respectively.
However the separation boundaries between adjacent classes have been shifted toward larger values: τ µ = 2 −20 (in-between microscopic and mesoscopic lengths) and τ m = 2 −6 (in-between mesoscopic and macroscopic lengths). The motivation stems from the need to start, over microscopic and mesoscopic length ranges, with computed integrand values over large enough sets of equally spaced quadrature knots. This enables a preliminary analysis of the complexity of the integral to be solved similar to that detailed below in the case of macroscopic ranges.
Over ranges of microscopic length, the minimal set of equally spaced quadrature knots is to consist of 2 × 2 + 1 = 5 elements. Over such a set, a two-term composite Simpson rule provides the reference quadrature sum output characterized by the highest possible algebraic degree of precision. The alternative knot grouping which consists of one central triplet for Simpson rule and two lateral doublets for composite trapezoidal rule yields a secondary composite quadrature sum output. From the difference of the two quadrature outputs and use of the triangle inequality rule, an associated error estimate results as a superposition of integrand curvatures at the three inner knots.
Over ranges of mesoscopic length, the use of a set of 4 × 2 + 1 = 9 equally spaced quadrature knots yields a reference quadrature sum from a two-term composite 4-interval Newton rule. Combination of a 4-interval Newton rule over the central knot quintuple with two lateral triplets for composite Simpson rule yields a secondary quadrature rule. Similar to the microscopic case, quadrature error estimates are obtained as a superposition of integrand curvatures at the seven inner knots.
Precision loss originating in cancellation by subtraction over microscopic and mesoscopic ranges. As noticed in [1], the specification of the integration domain ends, which enters the standard input formulation provided by the equation (1), results in an unavoidable precision loss due to the cancellation by subtraction involved in the computation of the integration domain length as the difference of two (arbitrary) machine numbers sharing a set of common most significant digits.
If, however, the integration domain length is specified to machine accuracy either separately or within the definition of the integral to be solved, then no precision loss originating in cancellation by subtraction will happen.

Enhanced accuracy Clenshaw-Curtis Quadrature
The recent monograph [10] provides a detailed rigorous presentation of the attractive mathematical properties of the CC quadrature. In the frame of the BAAQ, the selection of the CC quadrature sums for the approximation of the Riemann integrals (1) over macroscopic integration ranges is motivated by several features of the spanning Chebyshev polynomials: closeness to the polynomials of the best approximation, symmetry properties, easiness of the analytic computation of the quadrature weights for different weight functions [4,9,[11][12][13].
The interpolatory polynomial of the method is given by the truncated Chebyshev series expansion L φ n (y) equates φ(y) at the set of (n + 1) standard reduced CC quadrature knots {y n j = cos( jπ/n) | j = 0, 1, . . . , n}.
The computation of the Chebyshev expansion coefficients in (3) is the most computer intensive task of the method. The elucidation of the binary tree structures unveiling inheritance properties in the family of spanning Chebyshev polynomials of degrees n = 2 m (m = 1, 2, . . . ) [2] allowed the derivation of a fast algorithm for the computation of these coefficients inside which the number of multiplications is kept at the lowest possible value.
The scrutiny of the derived algorithm has shown that the CC quadrature spanned by the 32-nd degree Chebyshev polynomial provides the best CC solution (rich integrand sampling and still economical algorithm). An undesirable feature of the reduction procedure (4) to the standard integration domain [−1, 1] is the sizeable precision loss, due to the cancellation by subtraction, of the floating point values of the distances in-between neighbouring quadrature knots for the overwhelming fraction of the emerging pairs. This unfavourably impinges on the floating point accuracy of the divided differences entering the integrand profile (IP) analysis enabling Bayesian inferences [14].
The solution which avoids the above-mentioned precision loss is to use modified reduced abscissas (MRA), which are defined as distances of the standard reduced abscissas y n j to the nearest integration domain ends. The implementation of the MRA computation was done to machine accuracy. The computation of the distances between neighbouring MRA abscissas can be straightforwardly done, without precision loss by subtraction.
There is a set of 17 MRA entering the CC-32: The price to be paid for this accuracy increase is the separate computation and storage into two separate vectors of the integrand values at MRA over the left and the right halves of an integration

Flow chart of early Bayesian inference
There are two main results established in this subsection: (1) the early Bayesian inference enables consistent accommodation of the standard AAQ approach within the BAAQ and (2) robust termination criteria are defined which efficiently stop the computations while returning optimum output for the input integral (with a problem defined accuracy ceiling which supersedes an unreachable ceiling following from the user accuracy requirements). The flow chart below summarizes the eight steps of the proposed early Bayesian inference.
Step 6: GD diagnostic refinement based on the outputs got for the Chebyshev series expansion coefficients -GD is changed to EOC iff negligible highest label even-rank and odd-rank CC-32 coefficients -GD is changed to IC iff either (1) monotonicity is infringed for suitably chosen binary tree structure dependent subsets of CC-32 coefficients, or (2) non-convergence or slow convergence is established to occurr for the CC-32 coefficients -If (IC) then define offending IP abscissas and refine IC diagnostic at these abscissas; subdivide current subrange into diagnostic-dependent finer subranges elseif (GD) then proceed along the standard AAQ scheme.

Two illustrative examples
The accumulated empirical evidence plays a very important role in the decision taking process based on Bayesian inference. The data reported in [2] pointed to three different patterns of behaviour of the magnitudes of the Chebyshev expansion coefficients: fast convergence (suggesting rapid end of computations), moderate convergence (suggesting output improvement based on subrange bisection along the lines of the standard AAQ scheme), and ill-conditioning (irregular behaviour inside both the even-rank and the odd-rank coefficient subsets which is asking for IP analysis to resolve the offending integrand features).
The above-mentioned framework is refined by the scrutiny of the dependence of the even-rank (e) and odd-rank (o) Chebyshev expansion coefficients on the position of a singularity with respect to the ends of the current integration domain (figure 2).  In the left figure, the behaviour of the "e" and "o" coefficient subsets over the range [0, 1] is irregular due to the presence of an inner singularity at x s , in agreement with [2]. In this case, an ill-conditioning diagnostic is issued at step 8 of Sec. 4.2, which is asking for the localization of the offending integrand feature to machine accuracy. However, for the ranges [x s , 1] and [x s , √ 3/2], both characterized by endpoint singularities, a regular variation of both the "e" and "o" coefficient subsets is noticed and this asks in the decision path for the activation of a convergence acceleration procedure.
The data in figure 2-right are given in semi-logarithmic scale. The behavioural patterns over the inner derivative singularity range [−1, 1] as well as over the outer derivative singularity ranges [0, 1] and [−1, −0.75] are consistent with the data reported in Fig. 2-left and in [2]. The irregular behaviour pattern over the range [−1, 1] asks for localization of the offending abscissa to machine accuracy. Over the ranges [−1, x d ] and [−0.75, x d ], characterized by endpoint derivative singularities, the noticed slowly converging monotonic behaviours of the coefficient subsets ask for the activation of the convergence acceleration procedure. Far from singularity (the other two cases), the good convergence properties shift the decision path to the standard AAQ scheme.
In both the left and right figures 2, the smaller the extension of the integration domain with an endpoint singularity, the larger are the ranges of variation of the two Chebyshev expansion coefficient subsets.

Conclusions
The present report discusses a Bayesian automatic adaptive quadrature (BAAQ) solution for numerical integration which is simultaneously robust, reliable, and efficient, yielding maximum possible output accuracy in numerical experiments under arbitrary behaviour of the integrand function.
An essential ingredient of the solution is the multiscale approach.