Hybrid Symbolic-Numeric Integration in Multiple Dimensions via Tensor-Product Series 

Orlando A. Carvajal
Frederick W. Chapman
Keith O. Geddes 

Symbolic Computation Group, School of Computer Science
University of Waterloo, Waterloo, ON, N2L 3G1, Canada 

Abstract 

We present a new hybrid symbolic-numeric method for the fast and accurate evaluation of definite integrals in multiple dimensions. This method is well-suited for two classes of problems: (1) analytic integrands over general regions in two dimensions, and (2) families of analytic integrands with special algebraic structure over hyperrectangular regions in higher dimensions. 

 

The algebraic theory of multivariate interpolation via natural tensor-product series was developed in the doctoral thesis by Chapman, who named this broad new scheme of bilinear series expansions ”Geddes series” in honour of his thesis supervisor. This paper describes an efficient adaptive algorithm for generating bilinear series of Geddes-Newton type and explores applications of this algorithm to multiple integration. We will present test results demonstrating that our new adaptive integration algorithm is effective both in high dimensions and with high accuracy. For example, our Maple implementation of the algorithm has successfully computed nontrivial integrals with hundreds of dimensions to 10-digit accuracy, each in under 3 minutes on a desktop computer. 

 

Current numerical multiple integration methods either become very slow or yield only low accuracy in high dimensions, due to the necessity to sample the integrand at a very large number of points. Our approach overcomes this difficulty by using a Geddes-Newton series with a modest number of terms to construct an accurate tensor-product approximation of the integrand. The partial separation of variables achieved in this way reduces the original integral to a manageable bilinear combination of integrals of essentially half the original dimension. We continue halving the dimensions recursively until obtaining one-dimensional integrals, which are then computed by standard numeric or symbolic techniques. 

Introduction 

The hybrid symbolic-numeric method presented in this paper is well-suited for two classes of integration problems: 

 

 

Summary of the Method 

  • The method is based on recursive dimension-halving via separation of variables.
 

  • Multidimensional integrals are reduced to one-dimensional integrals using Geddes-Newton series expansions to generate accurate tensor-product approximations of the integrand.
 

  • Test results show that our new adaptive integration algorithm is effective both with high accuracy and in high dimensions.
 

Tensor-Product Series 

Tensor Products 

  • A tensor product is a finite sum of terms where each term is a product of univariate functions:
 

s[n](x, y) = Sum(g[i](x)*h[i](y), i = 1 .. n) 

 

  • There are many examples in mathematics:
 

exp(x+y) = exp(x)*exp(y) 

cos(x+y) = cos(x)*cos(y)-sin(x)*sin(y) 

(x+y)^2 = x^2+2*x*y+y^2 

  • The rank of a tensor product is the minimum number of terms among all the equivalent representations.
 

  • It is a natural tensor product when the factors of each term can be derived from the original function by a finite linear combination of linear functionals.
 

 

The Splitting Operator 

  • For a bivariate function  fand a point a, bwhere  f(a, b) <> 0, the splitting operator Upsilon[a, b]is defined by
       Upsilon[a, b]*f(x, y) = f(x, b)*f(a, y)/f(a, b).
 

  • The point a, b is called a splitting point.
 

  • Note: Upsilon[a, b] splits  f  into a rank-one tensor product.
 

  • Upsilon[a, b]*f(x, y)  interpolates  f(x, y) on the lines x = a and y = b , which means that  Upsilon[a, b]*`≡`(f(a, y), f(a, y))  and  Upsilon[a, b]*`≡`(f(x, b), f(x, b)).
 

  • These properties allow us to generate a natural tensor-product series simply by iterating the splitting operator while varying the choice of splitting point.
 

 

Geddes-Newton Series Expansions 

  • Not every bivariate function  f  can be expressed as a tensor product of finite rank.
 

  • For an approximation in  R = Typesetting:-delayCrossProduct([a, b], [c, d]) where  f  is continuous, we can uniformly approximate  f  by a rank-ntensor product:
 

s[n](x, y) = Sum(g[i](x)*h[i](y), i = 1 .. n) 

  • Define the n^th remainder :   r[n] = f-s[n] .
 

  • The following simple algorithm generates a natural tensor-product series expansion  s[n]  of a given function  f  on R .
 

Algorithm NTPS 

  • . Set  r[0] := f  and initialize i := 1 .
 

  • . Whiledelta < LinearAlgebra:-Norm(r[i-1], infinity) , iterate the following:
 

  • ) Choose `in`(a[i], b[i], R) with  abs(r[i-1](a[i], b[i])) = LinearAlgebra:-Norm(r[i-1], infinity) .
 

  • ) Let  r[i] := r[i-1]-Upsilon[a[i], b[i]]*r[i-1]  and  i := i+1 .
 

  • . Let  n := i-1  and return  s[n] := f-r[n] to approximate  f  with uniform error LinearAlgebra:-Norm(r[n], infinity) <= delta  on R .
 

 

  • Refining terminology introduced by Chapman in his doctoral thesis, we call  s[n]  the Geddes-Newton series expansion of  f  to  n  terms with respect to the splitting points  ({a[i], b[i]})[i = 1]^n .
 

Example 1 

Consider the following function on  Typesetting:-delayCrossProduct([0, 1], [0, 1]) : 

f(x, y) = exp(x^2*y^2)*cos(x+y) 

The Geddes-Newton series to 3 terms is as follows, where the splitting points are  a, a  for  a = 1, 0, .616336 : 

s[3](x, y) = Sum(c[i]*g[i](x)*h[i](y), i = 1 .. 3) 

where 

c[1] = -.884017 

g[1](x) = exp(x^2)*cos(x+1) 

h[1](y) = exp(y^2)*cos(y+1) 

c[2] = .794868 

g[2](x) = cos(x)+.477636*exp(x^2)*cos(x+1) 

h[2](y) = cos(y)+.477636*exp(y^2)*cos(y+1) 

c[3] = -9.83284 

g[3](x) = exp(.379870*x^2)*cos(x+.616336)-.356576*exp(x^2)*cos(x+1)-.623342*cos(x)
g[3](x) = exp(.379870*x^2)*cos(x+.616336)-.356576*exp(x^2)*cos(x+1)-.623342*cos(x)
g[3](x) = exp(.379870*x^2)*cos(x+.616336)-.356576*exp(x^2)*cos(x+1)-.623342*cos(x) 

h[3](y) = exp(.379870*y^2)*cos(y+.616336)-.356576*exp(y^2)*cos(y+1)-.623342*cos(y)
h[3](y) = exp(.379870*y^2)*cos(y+.616336)-.356576*exp(y^2)*cos(y+1)-.623342*cos(y)
h[3](y) = exp(.379870*y^2)*cos(y+.616336)-.356576*exp(y^2)*cos(y+1)-.623342*cos(y) 

Observations 

  • The splitting points in Example 1 are all located on the diagonal  y = x of the unit square, as illustrated in the following figure.
 

Plot 

  • Note:  s[3](x, y)  agrees exactly with  f(x, y)  at all points on the horizontal and vertical lines through the three splitting points!
 

  • With this strong interpolating property, even with only three splitting points we find that the maximum error over the unit square satisfies:  `≈`(LinearAlgebra:-Norm(f-s[3], infinity), `*`(3.0, `^`(10, -3))) .
 

  • If we request a series accurate to 15 digits for the function  f(x, y)  defined in Example 1, we find that only 12 terms are required (i.e., 12 splitting points). In other words,  LinearAlgebra:-Norm(f-s[12], infinity) < `*`(5.0, `^`(10, -15)) .
 

  • The distribution of the 12 splitting points for this particular function is illustrated in the following figure.
 

Plot 

Convergence of the Series 

  • Note similarities with the case of univariate polynomial interpolation via the Newton interpolation series.
 

  • In both cases, addition of a new term in the series corresponds to the addition of a new interpolation point.
 

  • However, in the case of Geddes-Newton series, the interpolation property holds not only at the new point but at all points on two lines passing through that point.
 

  • Intuitively, for any function  f  which is continuous in the unit square U, if the splitting points are dense along the diagonal then it is reasonable to expect that  s[n]  might converge uniformly to  f  on U as proc (n) options operator, arrow; infinity end proc.
 

  • In other words, by choosing well distributed splitting points on the diagonal, we can cover the whole unit square with small grid cells.  Recall: The remainder r[n] vanishes on the boundary of each cell .
 

  • Formal convergence results will appear in a separate paper.
 

A Brief Glance at Proof Technique 

  • The easiest case is when  f  is analytic in a sufficiently large region.
 

  • The convergence theorem and its proof are similar to the classical case where  f(z) is a function of a single complex variable and s[n](z) denotes the interpolating polynomial.
 

  • Recall the univariate case (Ref:  P.J. Davis, Interpolation and Approximation), as summarized below.
 

  • Basic tool: A contour integral representation of  f(z)-s[n](z) .
 

  • Let  f  be analytic on  `⊂`(S, complex)  and suppose that the interpolation points satisfy  `in`(a[i], `⊂`(A, S)) .
 

  • Let  s[n](z)  be the polynomial interpolating  f(z) at points {a[i]}.
 

  • alpha := diam(A),  delta := dist(A, `∂`(S)).
 

  • Then for all  `in`(z, A) ,abs(f(z)-s[n](z)) <= M*(alpha/delta)^n where  M = Typesetting:-delayCrossProduct(abs(`∂`(S))*`/`(2*Pi*delta), LinearAlgebra:-Norm(f, `∂`(S)))  .
 

  • Conclusion:  If  alpha < delta  then
 

limit(LinearAlgebra:-Norm(f-s[n], A), n = infinity) = 0 . 

 

Sketch 

  • For our case, we consider a bivariate function  f(w, z)  of two complex variables and  s[n](w, z)  denotes the Geddes-Newton series expansion of  f  to  n  terms with respect to the splitting points  {a[i], b[i]} .
 

  • Notation:   `in`(z, T)  and  S, `in`(b[i], `⊂`(B, T)) .
 

  • We derive a contour integration formula for the remainder  r[n] = f-s[n] :  


 

  • Here,  p[n](w) = product(w-a[k], k = 1 .. n)  and  q[n](z) = product(z-b[k], k = 1 .. n) .
 

  • Then we can conclude that  LinearAlgebra:-Norm(r[n], Typesetting:-delayCrossProduct(A, B)) <= M*gamma^n*LinearAlgebra:-Norm(r[n], Typesetting:-delayCrossProduct(`∂`(S), `∂`(T)))  where the constant  M  depends on the geometry, and gamma depends on separation constants (like before).
 

  • If  gamma < 1 and if the  n^th boundary constant satisfiesLinearAlgebra:-Norm(r[n], Typesetting:-delayCrossProduct(`∂`(S), `∂`(T))) = o(gamma^(-n)) as proc (n) options operator, arrow; infinity end proc then  limit(LinearAlgebra:-Norm(f-s[n], Typesetting:-delayCrossProduct(A, B)), n = infinity) = 0 .
 

Integration in Two Dimensions 

 

 

 

 

 

 

 

Transformation to Canonical 2-D Integral 

  • An integral with non-constant limits can be converted to a new integral with constant limits via a change of variables.
 

  • Specifically,
 

 

 

using the change of variables: 

 

x = a*(1-s)+b*s ;    y = c(x(s))*(1-t)+d(x(s))*t  

 

  • Therefore limit attention to integrals on the unit square.
 

 

  • Making the integrand symmetric: For an integral on the unit square, we can always make the integrand symmetric, meaning  `≡`(f(x, y), f(y, x)).
 

 

  • Namely,  f = f[S]+f[A]  (symmetric part + anti-symmetric part) where
 

 

f[S](x, y) = `*`(f(x, y)+f(y, x), 1/2) ;     

 

  • Integral of  f[A] over the unit square is 0, which gives us:
 

Int(Int(f(x, y), y = 0 .. 1), x = 0 .. 1) = Int(Int(f[S](x, y), y = 0 .. 1), x = 0 .. 1) . 

 

  • Therefore we always have a symmetric integrand on the unit square.
 

Example 2 

Consider the double integral 

Int(Int(exp(x^2*y^2)*cos(x+y), y = 0 .. 1-x), x = 0 .. 1) 

Applying the change of variables to transform to the unit square, and then applying symmetrization, yields the new problem: 

Int(Int(F(s, t), t = 0 .. 1), s = 0 .. 1) 

where 

F(s, t) = 1/2*cos(-s-t+s*t)*(exp(s^2*(1-s)^2*t^2)*(1-s)+exp(s^2*(1-t)^2*t^2)*(1-t))
F(s, t) = 1/2*cos(-s-t+s*t)*(exp(s^2*(1-s)^2*t^2)*(1-s)+exp(s^2*(1-t)^2*t^2)*(1-t)) 

Result of computation: 

  • Geddes-Newton series for  F(s, t)was computed to 3 terms
 

  • The series was integrated (applying one-dimensional quadratures)
 

  • Computed result:  0.385433
 

  • Result is correct to 5 significant digits (with only 3 terms)
 

Linear Algebra Formulation 

  • An efficient implementation is achieved by exploiting symmetry.
 

  • The series resulting from our algorithm has the algebraic form:
 

s[n](x, y) = Sum(c[i]*(Sum(k[i, j]*f(x, b[j]), j = 1 .. i))*(Sum(l[i, j]*f(a[j], y), j = 1 .. i)), i = 1 .. n) 

where  c[i], k[i, j], l[i, j] <> 0  are real-valued constants, and  a[i], b[i]  is the splitting point used in the i^th iteration. 

  • Note that  {a[i]} = {b[i]} .
 

  • Represent the series using matrices and vectors as
 

s[n](x, y) = V(x)^T*L^T*P*D*L*V(y) 

where 

  • V(x) is the column vector of dimension  n  with elements  f(x, a[i])
 

  • D  is an  Typesetting:-delayCrossProduct(n, n)  diagonal matrix with elements  c[i] = 1/r[i-1](a[i], b[i])
 

  • Pis a permutation matrix (identity matrix if splitting points all on diagonal)
 

  • L = [l[i, j]]  is an  Typesetting:-delayCrossProduct(n, n)  unit lower triangular matrix.
 

  • Benefits of the matrix-vector representation:
 

  • The cost of evaluating  s[n]  and  r[n]  is reduced from  O(n^2)  to  O(n)  evaluations of the original function  f .
 

  • The factors can be computed by exploiting efficient linear algebra operations for matrix-vector multiplication and inner product.
 

  • We only need to perform  n  one-dimensional integrations of cross-sections of the original function:   f(x, a[i])  for   .
 

  • The net effect:  Nearly all the processing time is spent evaluating the integrand, which is as desired.
 

Computational Results in 2-D 

The results in the following table are for the 2-D integral: 

Int(Int(sin(8*Pi*x*(1-x)*y*(1-y)*(x-y)^2), y = 0 .. 1), x = 0 .. 1) 

Tol 

Time 

Result 

Pts 

RelErr 

`*`(5, `^`(10, -10)) 

0.3 

0.069551393139 

15 

`*`(2.7, `^`(10, -13)) 

`*`(5, `^`(10, -15)) 

3.7 

0.06955139313890799 

21 

`*`(1.9, `^`(10, -19)) 

`*`(5, `^`(10, -20)) 

5.3 

0.0695513931389079901727 

26 

`*`(9.9, `^`(10, -24)) 

`*`(5, `^`(10, -25)) 

7.2 

0.069551393138907990172712980 

33 

`*`(4.5, `^`(10, -31)) 

`*`(5, `^`(10, -30)) 

14.2 

0.06955139313890799017271297964487 

38 

`*`(2.1, `^`(10, -35)) 

  • Tol = requested accuracy tolerance
 

  • Time = CPU time in seconds
 

  • Result = computed value of the integral (to the number of digits requested)
 

  • Pts = number of splitting points
 

  • RelErr = relative error, compared with an accurate reference value
 

Integration in High Dimensions 

 

 

 

 

 . 

 

 

 

The steps of DART 

  • . Deconstruction:  Deconstruct F to obtain a bivariate function f(s, t) = u[0](a*(s+t)+b)  or   f(s, t) = u[0]((as+b)*(at+b)) .
 

  • . Approximation:  Approximate the symmetric function  f  by a rank-n tensor-product series s[n] in two variables on the unit square ([0, 1])^2 .
 

  • . Reconstruction:  By a change of variables, reconstruct the original function F from f  and produce a d-variable tensor-product approximation S[n] of F on([0, 1])^d from the bivariate approximation s[n] of  f on ([0, 1])^2 .
 

Computational Results in High Dimensions 

The results below are for the following families of d-dimensional integrands: 

 

 

F[17] = ln(1+(sum(x[i], i = 1 .. d))) 

F[6] = (sum(b[i]*x[i], i = 1 .. d))^3 

F[18] = 1/(1+(sum(2*x[i], i = 1 .. d))) 

F[7] = cos(sum(x[i], i = 1 .. d)) 

F[19] = sin(1/4*Pi*(sum(x[i], i = 1 .. d)))^2 

F[16] = exp(product(x[i], i = 1 .. d)) 

F[22] = (sum(x[i], i = 1 .. d))^2*cos(sum(x[i], i = 1 .. d))
F[22] = (sum(x[i], i = 1 .. d))^2*cos(sum(x[i], i = 1 .. d)) 

 

Note:  F[5] and F[6] contain coefficients b[i] which were assigned non-zero integer values in the range   

 

Here, the region of integration is the unit d-cube:  `in`(x[i], [0, 1]) for all i , and the requested accuracy tolerance is fixed at  `*`(5, `^`(10, -10)) . 

Fcn 

Dim 

Time 

Result 

RelErr 

 

 

 

 

`*`(1.9, `^`(10, -13)) 

F[6] 

16 

21.3 

-`*`(1.9250000000, `^`(10, 2)) 

`*`(1.6, `^`(10, -13)) 

F[7] 

512 

37.9 

-`*`(1.8045810938, `^`(10, -11)) 

`*`(3.1, `^`(10, -10)) 

F[16] 

256 

165.2 

`*`(1.0000000000, `^`(10, 0)) 

`*`(2.0, `^`(10, -14)) 

F[17] 

64 

147.0 

`*`(3.4940406596, `^`(10, 0)) 

`*`(2.8, `^`(10, -12)) 

F[18] 

16 

30.4 

`*`(5.9973550251, `^`(10, -2)) 

`*`(1.5, `^`(10, -11)) 

F[19] 

128 

99.5 

`*`(4.9999927298, `^`(10, -1)) 

`*`(6.4, `^`(10, -13)) 

F[22] 

32 

92.6 

-`*`(5.6249710526, 10) 

`*`(9.8, `^`(10, -12)) 

  • Fcn = the integrand function
 

  • Dim = the dimension d for that integrand;  d was chosen based on the criterion that for the given integrand, dimension 2*d would require > 3 min
 

  • Time = CPU time in seconds
 

  • Result = computed value of the integral (to the number of digits requested)
 

  • RelErr = relative error, compared with an accurate reference value
 

References 

O.A. Carvajal, A New Hybrid Symbolic-Numeric Method for Multiple Integration Based on Tensor-Product Series. Master's thesis, University of Waterloo, Waterloo, Ontario, Canada, 2004. 

 

O.A. Carvajal, F.W. Chapman and K.O. Geddes, Hybrid symbolic-numeric integration in multiple dimensions via tensor-product series. Proceedings of ISSAC'05, Manuel Kauers (ed.), ACM Press, New York, 2005, pp. 84-91. 

 

F.W. Chapman, Generalized Orthogonal Series for Natural Tensor Product Interpolation. Ph.D thesis, University of Waterloo, Waterloo, Ontario, Canada, 2003. 

 

P.J. Davis, Interpolation and Approximation. Blaisdell, New York, 1963. 

 

D.H. Thomas, A natural tensor product interpolation formula and the pseudoinverse of a matrix. Linear Algebra and Its Applications, 13 (3), 1976, pp. 239-250.