Optimal Contours for High-Order Derivatives
OOPTIMAL CONTOURS FOR HIGH-ORDER DERIVATIVES
FOLKMAR BORNEMANN AND GEORG WECHSLBERGERA bstract . As a model of more general contour integration problems we considerthe numerical calculation of high-order derivatives of holomorphic functions usingCauchy’s integral formula. Bornemann ( ) showed that the condition numberof the Cauchy integral strongly depends on the chosen contour and solved theproblem of minimizing the condition number for circular contours. In this paperwe minimize the condition number within the class of grid paths of step size h using Provan’s algorithm for finding a shortest enclosing walk in weighted graphsembedded in the plane. Numerical examples show that optimal grid paths yieldsmall condition numbers even in those cases where circular contours are knownto be of limited use, such as for functions with branch-cut singularities. . I ntroduction To escape from the ill-conditioning of difference schemes for the numericalcalculation of high-order derivatives, numerical quadrature applied to Cauchy’sintegral formula has on various occasions been suggested as a remedy (for a surveyof the literature, see Bornemann ). To be specific, we consider a function f that is holomorphic on a complex domain D (cid:51)
0; Cauchy’s formula gives f ( n ) ( ) = n !2 π i (cid:90) Γ z − n − f ( z ) dz ( )for each cycle Γ ⊂ D that has winding number ind ( Γ ; 0 ) =
1. If Γ is not carefullychosen, however, the integrand tends to oscillate at a frequency of order O ( n − ) with very large amplitude (Bornemann , Fig. ). Hence, in general, there ismuch cancelation in the evaluation of the integral and ill-conditioning returnsthrough the backdoor. The condition number of the integral is (Deuflhard andHohmann , Lemma . ) κ ( Γ , n ) = (cid:82) Γ | z | − n − | f ( z ) | d | z | (cid:12)(cid:12) (cid:82) Γ z − n − f ( z ) dz (cid:12)(cid:12) and Γ should be chosen as to make this number as small as possible. Equivalently,since the denominator is, by Cauchy’s theorem, independent of Γ , we have tominimize d ( Γ ) = (cid:90) Γ | z | − n − | f ( z ) | d | z | . ( ) Mathematics Subject Classification. E , D ; R , C . Without loss of generality we evaluate derivatives at z = Given an accurate and stable (i.e., with positive weights) quadrature method such as Gauss–Legendre or Clenshaw–Curtis, this condition number actually yields, by ≈ log κ ( Γ , n ) ,an estimate of the error caused by round-off in the last significant digit of the data (i.e., the function f ). a r X i v : . [ m a t h . NA ] J u l FOLKMAR BORNEMANN AND GEORG WECHSLBERGER F igure 1 . Path Γ with ind ( Γ ; 0 ) = h . Bornemann ( ) considered circular contours of radius r ; he found that there isa unique r ∗ = r ( n ) solving the minimization problem and that there are differentscenarios for the corresponding condition number κ ∗ ( n ) as n → ∞ : • κ ∗ ( n ) → ∞ , if f is in the Hardy space H ; • lim sup n → ∞ κ ∗ ( n ) (cid:54) M , if f is an entire function of completely regulargrowth which satisfies a non-resonance condition of the zeros and whosePhragmén–Lindelöf indicator possesses M maxima (a small integer).Hence, though those (and similar) results basically solve the problem of choosingproper contours for entire functions, much better contours have to be found for theclass H . Moreover, the restriction to circles lacks any algorithmic flavor that wouldpoint to more general problems depending on the choice of contours, such as thenumerical solution of highly-oscillatory Riemann–Hilbert problems (Olver ). In this paper, we solve the contour optimization problem within the moregeneral class of grid paths of step size h (see Fig. ; we allow diagonals to beincluded) as they are known from Artin’s proof of the general, homological versionof Cauchy’s integral theorem (Lang , IV. ). Such paths are composed fromhorizontal, vertical and diagonal edges taken from a (bounded) grid Ω h ⊂ D ofstep size h . Now, the weight function ( ), being additive on the abelian group ofpath chains, turns the grid Ω h into an edge-weighted graph such that each optimalgrid path W ∗ becomes a shortest enclosing walk (SEW); “enclosing” because we haveto match the winding number condition ind ( W ∗ ; 0 ) =
1. An efficient solution ofthe SEW problem for embedded graphs was found by Provan ( ) and serves asa starting point for our work.
Outline of the Paper.
In Section we discuss general embedded graphs in whichan optimal contour is to be searched for; we discuss the problem of finding ashortest enclosing walk and recall Provan’s algorithm. In Section we discuss someimplementation details and tweaks for the problem at hand. Finally, in Section we give some numerical examples; these can easily be constructed in a way thatthe new algorithm outperforms, by orders of magnitude, the optimal circles ofBornemann ( ) with respect to accuracy and the direct symbolic differentiationwith respect to efficiency. Taking the contour optimization developed in this paper as a model, Wechslberger ( ) hasrecently addressed the deformation of Riemann–Hilbert problems from an algorithmic point of view.
PTIMAL CONTOURS FOR HIGH-ORDER DERIVATIVES . C ontour G raphs and S hortest E nclosing W alks By generalizing the grid Ω h , we consider a finite graph G = ( V , E ) embeddedto D , that is, built from vertices V ⊂ D and edges E that are smooth curvesconnecting the vertices within the domain D . We write uv for the edge connectingthe vertices u and v ; by ( ), its weight is defined as d ( uv ) = (cid:90) uv | z | − n − | f ( z ) | d | z | . ( )A walk W in the graph G is a closed path built from a sequence of adjacent edges,written as (where ˙ + denotes joining of paths) W = v v ˙ + v v ˙ + · · · ˙ + v m v ;it is called enclosing the obstacle 0 if the winding number is ind ( W ; 0 ) =
1. The setof all possible enclosing walks is denoted by Π . As discussed in § , the conditionnumber is optimized by the shortest enclosing walk (not necessarily unique) W ∗ = argmin W ∈ Π d ( W ) where, with W = v v ˙ + v v ˙ + · · · ˙ + v m v and v m + = v , the total weight is d ( W ) = m ∑ j = d ( v j v j + ) .The problem of finding such a SEW was solved by Provan ( ): the idea is thatwith P u , v denoting a shortest path between u and v , any shortest enclosing walk W ∗ = w w ˙ + w w ˙ + · · · ˙ + w m w can be cast in the form (Provan , Thm. ) W ∗ = P w , w j ˙ + w j w j + ˙ + P w j + , w for at least one j . Hence, any shortest enclosing walk W ∗ is already specified byone of its vertices and one of its edges; therefore W ∗ ∈ ˜ Π = {P u , v ˙ + vw ˙ + P w , u : u ∈ V , vw ∈ E } .Provan’s algorithm finds W ∗ by, first, building the finite set ˜ Π ; second, by removingall walks from it that do not enclose z =
0; and third, by selecting a walk fromthe remaining candidates that has the lowest total weight. Using Fredman andTarjan’s ( ) implementation of Dijkstra’s algorithm to compute the shortestpaths P u , v , the run time of the algorithm is known to be (Provan , Corollary ) O ( | V | | E | + | V | log | V | ) . ( ) . I mplementation D etails We restrict ourselves to graphs Ω h given by finite square grids of step size h ,centered at z = D . Since Provan’s algorithm just requires an embedded graph but nota planar graph, we may add the diagonals of the grid cells as further edges tothe graph (see Fig. ). For such a graph Ω h , with or without diagonals, we have These diagonals increase the number of possible slopes which results, e.g., in improved approxi-mations of the direction of steepest descent at a saddle point of d ( z ) (Bornemann , § ) or in a fasterU-turn around the end of a branch-cut, see Fig. . The latter case leads to some significant reductionsof the condition number, see Fig. . FOLKMAR BORNEMANN AND GEORG WECHSLBERGER Ai ( z ) exp ( ( + z ) )( − z ) J ( z ) F igure 2 . W ∗ (red) vs. W v ∗ (blue): the color coding shows the size of log d ( z ) ; withred for large values and green for small values. The smallest level shown is thethreshold, below of which the edges of W ∗ do not contribute to the first couple ofsignificant digits of the total weight. The plots illustrate that W ∗ and W v ∗ differtypically just in a small region well below this threshold; consequently, both walksyield about the same condition number. On the right note the five-leaved cloverthat represents the combination of algebraic and essential singularity at z = − | V | = O ( h − ) and | E | = O ( h − ) so that the complexity bound ( ) simplifies to O ( h − log h − ) . . . Edge Weight Calculation.
Using the edge weights d ( uv ) on Ω h requires toapproximate the integral in ( ). Since not much accuracy is needed here, a simpletrapezoidal rule with two nodes is generally sufficient: d ( uv ) = (cid:90) uv | z | − ( n + ) | f ( z ) | d | z | = | u − v | ( d ( u ) + d ( v )) + O ( h ) = ˜ d ( uv ) + O ( h ) with the vertex weight d ( z ) = | z | − ( n + ) | f ( z ) | . ( )Although ˜ d ( uv ) will typically have an accuracy of not more than just a few bitsfor the rather coarse grids Ω h we work with, we have not encountered a singlecase in which a more accurate computation of the weights would have resulted ina different SEW W ∗ . . . Reducing the size of ˜ Π . As described in Section , Provan’s algorithm startsby building a walk for every pair ( v , e ) ∈ V × E and then proceeds by selectingthe best enclosing one. A simple heuristic, which worked well for all our test cases,helps to considerably reduce the number of walks to be processed: Let v ∗ = argmin v ∈ V d ( v ) Recall that optimizing the condition number is just a question of order of magnitude but not ofprecise numbers. Once the contour Γ has been fixed, a much more accurate quadrature rule will beemployed to calculate the integral ( ) itself, see § . . PTIMAL CONTOURS FOR HIGH-ORDER DERIVATIVES and define W v ∗ as a SEW subject to the constraint W v ∗ ∈ ˜ Π v ∗ = {P v ∗ , u ˙ + uw ˙ + P w , v ∗ : uw ∈ E } .Obviously W ∗ and W v ∗ do not need to agree in general, as v ∗ does not have tobe traversed by W ∗ . However, since v ∗ is the vertex with lowest weight, bothwalks differ mainly in a region that has no, or very minor, influence on the totalweight and, consequently, also no significant influence on the condition number.Actually, W ∗ and W v ∗ yielded precisely the same total weight for all functions thatwe have studied (Fig. compares W ∗ and W v ∗ for two typical examples). Usingthat heuristic, the run time of Provan’s algorithm improves to O ( | E | + | V | log | V | ) because its main part reduces to applying Dijkstra’s shortest path algorithm justonce. In the case of the grid Ω h this bound simplifies to O ( h − log h − ) . . . Size of the Grid Domain.
The side length l of the square domain supporting Ω h has to be chosen large enough to contain a SEW that would approximate anoptimal general integration contour. E.g., if f is entire, we choose l large enoughfor this square domain to cover the optimal circular contour: l > r ∗ , where r ∗ is the optimal radius given in Bornemann ( ); a particularly simple choice is l = r ∗ . In other cases we employ a simple search for a suitable value of l bycalculating W ∗ for increasing values of l until d ( W ∗ ) does not decrease substantiallyanymore. During this search the grid will be just rescaled, that is, each grid usesa fixed number of vertices; this way only the number of search steps enters as anadditional factor in the complexity bound. . . Multilevel Refinement of the SEW.
Choosing a proper value of h is notstraightforward since we would like to balance a good approximation of a generallyoptimal integration contour with a reasonable amount of computing time. Inprinciple, we would construct a sequence of SEWs for smaller and smaller valuesof h until the total weight of W ∗ does not substantially decrease anymore. To avoidan undue amount of computational work, we do not refine the grid everywherebut use an adaptive refinement by confining it to a tubular neighborhood of thecurrently given SEW W ∗ (see Fig. ): : calculate W ∗ within an initial grid; : subdivide each rectangle adjacent to W ∗ into rectangles; : remove all other rectangles; : calculate W ∗ in the newly created graph.As long as the total weight of W ∗ decreases substantially, steps to are repeated.It is even possible to tweak that process further by not subdividing rectangles thatjust contain vertices or edges of W ∗ having weights below a certain threshold. Bygeometric summation, the complexity of the resulting algorithm is O ( H − log H − ) + O ( h − log h − ) where H denotes the step size of the coarsest grid and h = H /2 k the step sizeafter k loops of adaptive refinement. An analogous approach to the constrained W v ∗ -variant of the SEW algorithm given in § . reduces the complexity further to O ( H − log H − ) + O ( h − log h − ) , FOLKMAR BORNEMANN AND GEORG WECHSLBERGER FOLKMAR BORNEMANN AND GEORG WECHSLBERGER step step step step comparison comparisoninitial graphinitial walkfirst refined gridfirst refined walk F igure 3 . Multilevel refinement of W ⇤ ( f ( z ) = G ( z ) , n = ) would need: Fig. shows an example with the order n =
300 of differentiationbut accurate solutions using just about 200 nodes which is well below what thesampling condition would require for circular contours (Bornemann , § . ).Of course, trapezoidal sums would also benefit from some recursive device thathelps to neglect those nodes which do not contribute to the numerical result. . N umerical R esults Table displays condition numbers of SEWs W ⇤ on rectangular grids as com-pared to the optimal circles C r ⇤ for a couple of functions; Fig. shows some of thecorresponding contours. For entire f we observe that W ⇤ , like the optimal circle C r ⇤ ,automatically traverses the saddle points of d ( z ) . It was shown in Bornemann F igure 3 . Multilevel refinement of W ∗ ( f ( z ) = Γ ( z ) , n = ) which is close to the best possible bound O ( h − ) given by the work that would beneeded to just list the SEW. . . Quadrature Rule for the Cauchy Integral.
Finally, after calculation of theSEW Γ = W ∗ , the Cauchy integral ( ) has to be evaluated by some accurate numerical quadrature. We decompose Γ into maximally straight line segments,each of which can be a collection of many edges. On each of those line segments weemploy Clenshaw–Curtis quadrature in Chebyshev–Lobatto points. Additionallywe neglect segments with a weight smaller than 10 − times the maximum weightof an edge of Γ , since such segments will not contribute to the result withinmachine precision. This way we not only get spectral accuracy but also, in manycases, less nodes as would be needed by the vanilla version of trapezoidal sums on acircular contour: Fig. shows an example with the order n =
300 of differentiationbut accurate solutions using just about 200 nodes which is well below what thesampling condition would require for circular contours (Bornemann , § . ).Of course, trapezoidal sums would also benefit from some recursive device thathelps to neglect those nodes which do not contribute to the numerical result. PTIMAL CONTOURS FOR HIGH-ORDER DERIVATIVES r e l a t i v ee rr o r f ( z ) = (1 − z ) / , n = 10without diagonalswith diagonals r e l a t i v ee rr o r f ( z ) = (1 − z ) / , n = 300without diagonalswith diagonals F igure 4 . Illustration of the spectral accuracy of piecewise Clenshaw–Curtis quad-rature on SEW contours for a function with a branch-cut singularity. For larger n ,we observe a significant improvement by adding diagonals to the grid. We get tomachine precision for n =
10 and loose about two digits for n = digits for n = digits for n = , Thm. . ).T able 1 . Condition numbers for some f ( z ) : r ∗ are the optimal radii given inBornemann ( ); W ∗ was calculated in all cases on a 51 × l = r ∗ (in the last two cases l was found as in § . ). For 1/ Γ ( z ) , the peculiar order ofdifferentiation n = , Table ). In the last example, differentiation is for z = √ f ( z ) n κ ( W ∗ , n ) κ ( C r ∗ , n ) e z
300 1.1 1.0Ai ( z )
300 1.3 1.21/ Γ ( z )
300 1.7 1.61/ Γ ( z ) · · ( − z )
10 1.4 5.0 · exp ( ( + z ) )( − z ) J ( z )
100 7.2 · · T able 2 . CPU times for the examples of Table . Here t W ∗ and t W v ∗ denote the timesto compute W ∗ and W v ∗ and t quad denotes the time to approximate the integral ( )on such a contour by quadrature. (There is no difference between W ∗ and W v ∗ from the point of quadrature, see Fig. .) In the last example, differentiation is for z = √
2. The timings for the grids of size 25 ×
25 and 51 ×
51 match nicely the O ( h − log h − ) complexity for W ∗ and the O ( h − log h − ) complexity for W v ∗ . f ( z ) n grid t W ∗ t W v ∗ t quad e z
300 51 ×
51 4.4 · s 1.5 s 0.3 sAi ( z )
300 25 ×
25 2.1 · s 0.5 s 1.7 sAi ( z )
300 51 ×
51 4.0 · s 2.1 s 2.1 s1/ Γ ( z )
300 25 ×
25 2.0 · s 0.5 s 1.5 s1/ Γ ( z )
300 51 ×
51 3.6 · s 2.4 s 1.3 s1/ Γ ( z ) ×
51 3.6 · s 2.3 s 3.1 s ( − z )
10 51 ×
51 1.4 · s 5.9 s 0.2 sexp ( ( + z ) )( − z ) J ( z )
100 51 ×
51 7.0 · s 3.5 s 0.3 s FOLKMAR BORNEMANN AND GEORG WECHSLBERGER Ai ( z ) exp ( ( + z ) )( − z ) J ( z ) F igure 5 . W v ∗ (blue: Ω h without diagonals, magenta: Ω h with diagonals) vs. C r ∗ (cyan) for some examples of Table : the color coding shows the size of log d ( z ) ;with red for large values and green for small values. The smallest level shownis the threshold, below of which the edges of W v ∗ do not contribute to the firstsignificant digits of the total weight. . N umerical R esults Table displays condition numbers of SEWs W ∗ as compared to the optimalcircles C r ∗ for five functions; Table gives the corresponding CPU times and Fig. shows some of the contours. ( All experiments were done using hardware arithmetic. )The purpose of these examples is twofold, namely to demonstrate that:( ) the SEW algorithm matches the quality of circular contours in cases wherethe latter are known to be optimal such as for entire functions;( ) the SEW algorithm is significantly better than the circular contours in caseswhere the latter are known to have severe difficulties.Thus, the SEW algorithm is a flexible automatic tool that covers various classes ofholomorphic functions in a completely algorithmic fashion; in particular there isno deep theory needed to just let the computation run.In the examples of entire f we observe that W ∗ and W v ∗ , like the optimal cir-cle C r ∗ would do, traverses the saddle points of d ( z ) . It was shown in Bornemann( , Thm. . ) that, for such f , the major contribution of the condition numbercomes from these saddle points and that circles are (asymptotically, as n → ∞ )paths of steepest decent. Since W ∗ can cross a saddle point only in a horizontal,vertical, or (if enabled) diagonal direction, somewhat larger condition numbershave to be expected. However, the order of magnitude of the condition numberof C r ∗ is precisely matched. This match holds in cases where circles give a con-dition number of approximately 1, as well as in cases with exceptionally largecondition numbers, such as for f ( z ) = Γ ( z ) in the peculiar case of the order ofdifferentiation n = , § . ).For non-entire f , however, optimized circles will be far from optimal in general:Bornemann ( , Thm. . ) shows that the optimized circle C r ∗ for functions f from the Hardy space H with boundary values in C k , α yields a lower conditionnumber bound of the form κ ( C r ∗ , n ) (cid:62) cn k + α ; PTIMAL CONTOURS FOR HIGH-ORDER DERIVATIVES n c o nd i t i o nnu m b e r f ( n ) ( z ) = exp ( − / (1 + 8 z ))(1 − z ) / J ( z ), z = 1 / √ W ∗ C . n / F igure 6 . An example with essential and algebraic singularities: the conditionnumber of the Cauchy integral for exp ( ( + z ) )( − z ) J ( z ) for varyingorder n of differentiation at z = √
2; blue: optimal contour W ∗ in a 51 ×
51 gridgraph; green: circular contour with near optimal radius r = ≈ − √
2; red:prediction of the growth rate from Bornemann ( , Thm. . ). for instance, f ( z ) = ( − z ) gives κ ( C r ∗ , n ) ∼ · n . On the otherhand, W ∗ gives condition numbers that are orders of magnitude better than thoseof C r ∗ by automatically following the branch cut at ( ∞ ) .The latter example can easily be cooked-up to outperform symbolic differ-entiation as well: using Mathematica , the calculation of the n -th derivative of f ( z ) = exp ( ( + z ) )( − z ) J ( z ) at z = √ n =
23 but had to be stopped after more than a week for n = z = −
1, the W v ∗ version of the SEW calculates this n = digits in less than 4 s; whereas optimized circularcontours would give only about correct digits here (see Fig. ).While many more such numerical experiments would demonstrate that reason-ably small condition numbers are obtainable in general, the study of rigorouscondition number bounds for the SEW has to be postponed to future work.R eferences Bornemann, F.: , Accuracy and stability of computing high-order derivatives of analytic functionsby Cauchy integrals,
Found. Comput. Math. , – .Deuflhard, P. and Hohmann, A.: , Numerical analysis in modern scientific computing , second edn,Springer-Verlag, New York.Fredman, M. L. and Tarjan, R. E.: , Fibonacci heaps and their uses in improved network optimizationalgorithms,
J. Assoc. Comput. Mach. , – .Lang, S.: , Complex analysis , fourth edn, Springer-Verlag, New York.Olver, S.: , Numerical solution of Riemann–Hilbert problems: Painlevé II,
Found. Comput. Math. , – .Provan, J. S.: , Shortest enclosing walks and cycles in embedded graphs, Inform. Process. Lett. , – .Wechslberger, G.: , Automatic deformation of Riemann-Hilbert problems. arXiv:1206.2446 .Z entrum M athematik – M , T echnische U niversität M ünchen , M ünchen , G ermany E-mail address : [email protected]; [email protected] The software is provided as a supplement to the e-print version of this paper: arXiv:1107.0498arXiv:1107.0498