Quantifying the ill-conditioning of analytic continuation

Analytic continuation is ill-posed, but becomes merely ill-conditioned (although with an infinite condition number) if it is known that the function in question is bounded in a given region of the complex plane. In an annulus, the Hadamard three-circles theorem implies that the ill-conditioning is not too severe, and we show how this explains the effectiveness of Chebfun and related numerical methods in evaluating analytic functions off the interval of definition. By contrast, we show that analytic continuation is far more ill-conditioned in a strip or a channel, with exponential loss of digits of accuracy at the rate $\exp(-\pi x/2)$ as one moves along. The classical Weierstrass chain-of-disks method loses digits at the faster rate $\exp(-e\kern .3pt x)$.

assumptions being made explicit-and this is understandable, for in applications, often one has a sense of certain features of one's function without being able to pin them down precisely. In this paper, however, we wish to be completely explicit and show how certain natural regularizing assumptions lead to upper and lower bounds on the accuracy of analytic continuation.
Our regularizing assumption will be that f is not only analytic in Ω , but bounded. For example, we can take the bound to be 1 2 and consider the set of functions that are analytic in Ω and satisfy f Ω ≤ 1 2 . (The symbol · A always denotes the supremum norm over the set A.) If f ,f are two such functions, then f − f Ω ≤ 1. We shall show (Theorem 5.2) that if in addition f − f E ≤ ε, then for each z ∈ Ω \E, for some α(z) ∈ (0, 1) that depends on Ω , E, and z but not on f andf or ε. Another way to say the same thing is We may interpret (1.2) as follows: if we know a function f satisfying f Ω ≤ 1 2 to d digits on E, then it is determined to α(z)d digits at z. Even though our theorems are, of course, mathematical rather than computational results, we shall use the terminology of digits a good deal in discussing them, for this is an easy way to talk about logarithmic quantities.
This general framework may sound rather abstract. It becomes concrete when we consider the dependence of α(z) on z for particular choices of Ω and E, and after stating a basic lemma in Section 2, we shall focus on two choices that are particularly fundamental. The first is radial geometry, with analytic continuation outward from the unit disk E into a disk Ω of radius R > 1 (Section 3). In this setting analytic continuation is reasonably well-conditioned, with digits of accuracy being lost only linearly as |z| increases. This observation possibly goes back to Hadamard himself, and its numerical implications have been considered by various authors including Miller [20] and Franklin [12]. An intuitive way to understand the effect is to note that in the limiting case R → ∞, Liouville's theorem implies that f must be constant, so if we know f to accuracy ε on E, we know it to the same accuracy everywhere. The result for finite R (Theorem 3.1) can be derived from the Hadamard three-circles theorem [18]. The essence of the matter is that analyticity and boundedness at a large radius imply rapid exponential decrease of Taylor coefficients, hence good behavior at smaller radii.
Analytic continuation is much more difficult in the other geometry we focus on, which is linear (Section 4). Here we take Ω to be an infinite half-strip of half-width 1 (without loss of generality), and E as the end segment of the half-strip. As z moves away from E along the centerline of the strip, digits are lost exponentially as a function of distance, and we prove this by reducing the problem to the configuration of Lemma 2.1. At a point z that is 2π units away from the end, for example, the number of accurate digits has shrunk by a factor (π/4) exp(π 2 ) ≈ 15,000, so if you want to have 3 digits of accuracy at such a point, you'll need to start with 45,000 digits.
Our formulations are conformally invariant, and thus different regions Ω = C can be transplanted from one to another. In particular, analytic continuation in a disk and a half-strip are essentially equivalent problems, and the reason the half-strip is exponentially more difficult than the disk is that the conformal map that relates them is an exponential. Section 5 explores results for general regions (not necessarily simply connected) that follow from these observations, presenting theorems establishing the behavior asserted in the opening paragraphs of this introduction.
Along the way, we shall relate our results to numerical algorithms. Section 3 presents a simple method for numerical analytic continuation in a disk that approximately achieves the bounds indicated in Theorem 3.1, based on Taylor series on the disk, and in Section 6 we show that this method is implicit in Chebfun [9,26]. An algorithm of this kind was proposed by Franklin [12], and there is recent related work by Demanet and coauthors [4,7], among others. Further numerical algorithms and associated mathematical estimates for analytic continuation can be found in [6,8,13,14,15,16,17,20,22,23,25,27]. More generally, there is a large literature of numerical methods for ill-posed problems, which are often defined by partial differential or integral equations. One paper that speaks of the connection between analytic continuation and more general ill-posed problems defined by PDEs is [20].
In Section 7 we turn to the most famous algorithm of analytic continuation, which goes back to Weierstrass: marching Taylor expansions from one overlapping disk to another in a chain. We show that this method, if carried out numerically in the half-strip with a certain optimal choice of parameters, suffers exponential loss of accuracy at a rate 2e/π times faster than the optimal rate in a half-strip, so that if one marches 2π units down the half-strip, the number of accurate digits is divided by exp(2πe) ≈ 26,000,000.
Before turning to the details, we comment on the relationship between approximation methods based on multiple derivative values at a single point, such as chainof-disks continuation of Taylor series or Padé approximation [3], and methods based just on function values but at multiple points. Our formulations are of the latter form, but the two contexts are close. Thanks to the standard lemma of complex analysis known as Cauchy's estimate, knowing a function f on the unit disk to accuracy ε is approximately the same as knowing its Taylor coefficients c k to accuracy ε, and more generally, if f is known to accuracy ε on the closed disk of radius r, this is approximately the same as knowing its Taylor coefficients c k to accuracy ε r −k .

A lemma
Our results are based on the following lemma, the Hadamard three-lines theorem, a more general form of which can be found in [24,Thm. 12.8]. Numerical algorithms for this geometry are discussed in [14].
Proof We begin by noting that without loss of generality, we may suppose that h is analytic on the closed set S. If not, we could restrict attention to a smaller domain δ ≤ Re w ≤ 1 − δ for δ < 0 and then take the limit δ → 0. Set ν = − log ε, implying e −ν = ε. The function e −ν w h(w) is analytic in S and bounded in absolute value by ε for Re w = 0 and also for Re w = 1. Therefore, by the maximum modulus principle as qualified in the next paragraph, it is bounded by ε for all w ∈ S. Thus |h(w)| ≤ ε e ν Re w = e (1−Re w) log ε , as required.
The qualification just mentioned is that the maximum modulus principle does not apply to arbitrary functions on an unbounded domain with a gap in the boundary at ∞. However, this function is known to be bounded in an infinite strip, and in such a situation, according to a Phragmén-Lindelöf theorem [18,Thm. 18.1.4], the maximum modulus principle applies after all.
For the converse, it is enough to consider the function h(w) = ε e νw . ⊓ ⊔ .1 shows our first fundamental geometry, one that has been considered by a number of authors. We take E to be the closed unit disk and Ω as the open disk of radius R > 1. As always, our concern is obtaining bounds on |g(z)|, z ∈ Ω \E, for an analytic function g (=f − f ) satisfying g Ω ≤ 1 and g E ≤ ε. Here is the result, essentially the Hadamard three-circles theorem, with · r denoting the supremum norm over {z : |z| < r}.

Analytic continuation in a disk
Theorem 3.1 Given R > 1, let g be analytic in Ω = {z ∈ C : |z| < R} with g R ≤ 1 and g 1 ≤ ε ∈ (0, 1). Then for any z with 1 < |z| < R, Conversely, for an infinite sequence of values ε converging to 0, there are functions g satisfying the given conditions for which the inequalities (3.1) hold as equalities for all z with 1 < |z| < R.
Proof As in the proof of Lemma 2.1, we begin by noting that without loss of generality, we may suppose that g is analytic in the closed domain |z| ≤ R. We transplant the problem to the infinite strip S of the last section by defining h(w) = g(z) for z = R w = e w log R , hence w = log z/ log R (it doesn't matter which branch of log z is used, as they all lead to the same estimate). By the lemma, we have log For the converse result, consider g(z) = (z/R) n with n = − log ε/ log R. For an infinite sequence of values ε → 0, n is an integer, and in these cases g is an analytic function with the required properties.
⊓ ⊔ If f is known to d digits on E, the number of digits determined in Ω falls off smoothly to 0 at the outer boundary.
In words, we can describe Theorem 3.1 as follows. If f is known to d digits for |z| ≤ 1, the number of digits determined for |z| = r diminishes to 0 as r → R; as a function of log r, the loss of digits is linear.
It is easy to outline an algorithm for analytic continuation, based on a finite Taylor series of length n ≈ − log ε/ log R, that achieves approximately the accuracy promised in Theorem 3.1. From approximate values for |z| = 1 with error at most ε of a function f with f Ω ≤ 1 2 , we compute approximationsc k ≈ c k to the Taylor coefficients {c k } of f for 0 ≤ k ≤ n with |c k − c k | ≤ ε. That this is possible follows from Cauchy's estimate applied on the circle |z| = 1; in practice, we sample f on a grid of N ≫ n roots of unity and use the Fast Fourier Transform [2]. By Cauchy's estimate applied now for |z| = R, the Taylor coefficients of f satisfy |c k | ≤ 1 2 R −k . If we definẽ Our choice n ≈ − log ε/ log R implies ε ≈ R −n , and thus these two sums are both of size on the order of (|z|/R) n . We therefore achieve, as required, with the equality in the middle holding since the definition (3.2) implies R α(z) = R/|z|. In the algorithm just described, the regularization occurred when we took the series (3.3) to be finite rather than infinite. In this geometry, the ill-posedness of analytic continuation resides in the fact that as k → ∞, powers z k have unbounded discrepancies of absolute value between one radius |z| and another. For a fascinating analysis of the implications of such behavior in the context of computation of Taylor coefficients, see [5]. The importance of truncating a series for numerical analytic continuation was recognized at least as early as [19], and a detailed error analysis of an algorithm with this flavor can be found in [12].
We have described the algorithm as a process for working with a function known to accuracy ε on the unit disk. One way to obtain such data is to take Taylor polynomials of f of successively higher degrees (in exact arithmetic), in which case (3.1) amounts to the statement that the convergence of Taylor polynomials to f (z) for |z| < R is exponential at a rate of order (|z|/R) n .

Analytic continuation in a half-strip
Our second fundamental geometry is shown in Figure 4.1. Now E is the complex interval (−i, i) and Ω is the half-strip H of points z = x + iy with x > 0, −1 < y < 1. Our aim is to analytically continue a function f from the end segment of the strip to real positive values x.
To apply Lemma 2.1 in this geometry, we need to map H conformally to the infinite strip S of Section 2 in such a way that the corner points z = ±i map to the infinite vertices w = ±i∞ and z = ∞ maps to w = 1. We can construct such a map by composing three simpler maps. First, u = sinh(πz/2) maps H to the right halfplane with distinguished points ±i and ∞. Next, v = (i − u)/(i + u) maps the right half-plane to the upper half-plane with distinguished points 0, ∞, and −1. Finally, w = (−i/π) log v maps the upper half-plane to S as required. Combining these steps, we find that the map from the half-strip H in the z-plane to the infinite strip S in the w-plane is given by The properties of this map are such as to reveal that analytic continuation into a half-strip where a function is known to be bounded, while possible in principle, is so ill-conditioned as to be generally infeasible in practice. Instead of digits of accuracy being lost linearly, they are lost exponentially, as emphasized in Figure 4.1. To see how this comes about, consider the situation in which z is a real number x > 0. The leading terms in the asymptotics for u, v, and w for large x give us Since w is exponentially close to 1, Lemma 2.1 implies that the number of digits of accuracy will be multiplied by an exponentially small factor ∼ (4/π) exp(−πx/2). A more careful analysis sharpens (4.2) to from which we get the following theorem: Theorem 4.1 Let g be analytic in the half-strip H with g H ≤ 1 and g E ≤ ε for some ε ∈ (0, 1). Then for any x > 0,

5)
Conversely, for any ε ∈ (0, 1), there is a function g satisfying the given conditions for which, for all x ≥ 0, Proof This follows from Lemma 2.1 by the conformal transplantation (4.1), using the bounds (4.3). ⊓ ⊔

General geometries
For more general geometries than the disk or the half-strip, let us now justify the claims made in the introduction. First is the ill-posedness statement of the opening paragraph, which as usual we formulate as an assertion about an analytic function g =f − f .

Theorem 5.1
Let Ω be a connected open region of the complex plane C and let E be a bounded nonempty continuum in Ω whose closure E does not enclose any points of Ω \E. Let g be an analytic function in Ω ∪ E satisfying g E ≤ ε for some ε > 0. This condition implies no bounds whatsoever on the value of g at any point z ∈ Ω \E.
Proof Given z ∈ Ω \E and a complex number M, we shall show there is a polynomial p such that p E ≤ ε and p(z) = M. LetẼ denote the compact set consisting of E together with all points enclosed by this set; thus the complement ofẼ in the complex plane is connected. By assumption, z ∈Ẽ. According to Runge's theorem [24,Thm. 13.7], there is a polynomial q such that q Ẽ ≤ ε/2 and |q(z) − M| ≤ ε/2. Now define p(ζ ) = q(ζ ) + M − q(z). We readily verify p(z) = M and | p(ζ )| ≤ |q(ζ )| + |M − q(z)| ≤ ε for all ζ ∈ E.

⊓ ⊔
The well-posedness statement of the introduction is (1.1), which we formulate as follows.
Theorem 5.2 Let Ω , E, ε, and g be as in Theorem 5.1, but now with g additionally satisfying g Ω ≤ 1, and let z be a point in Ω \E. Assume that the boundary of E is piecewise smooth (a finite union of smooth Jordan arcs). Then there is a number α ∈ (0, 1), independent of g though not of z, such that for all ε > 0, |g(z)| ≤ ε α(z) . Theorem 5.2 asserts that analytic continuation in the presence of a boundedness condition is a well-posed problem in the sense that there is a unique solution depending continuously on the data, but this does not mean that its condition number is finite. A finite condition number would correspond to |g(z)| shrinking linearly with ε, that is, to a value α = 1 in (5.1), or more generally to the bound as ε → 0 for some constant C depending on z but not g. But Theorem 5.2 only gives α < 1, and in fact, we now show that α = 1 cannot occur except in trivial cases with Ω = C. can never hold unless Ω is the whole complex plane C.
Proof If Ω is all of C, we may be in the trivial situation mentioned in the introduction, where Liouville's theorem implies that g is constant. (This will be true, for example, if Ω consists of C with a finite set of points removed. It won't be true if Ω consists of C with some arcs removed.) On the other hand suppose there is a point z 0 ∈ C disjoint from Ω . Then there is a closed disk ∆ about z 0 that is disjoint from Ω . By the conformal map 1/(z − z 0 ), we may transplant the problem so that Ω and E are bounded. For any z ∈ Ω \E, as in the proof of Theorem 5.1, Runge's theorem ensures that there is a polynomial p such that p(z) = 0 and Re p(ζ ) ≤ −1 for all ζ ∈ E. Choose M > 0 such that Re p(ζ ) ≤ M for all ζ ∈ Ω . Then q(z) = p(z) − M has real parts ≤ −(1 + M) for ζ ∈ E, −M at z, and ≤ 0 for ζ ∈ Ω . Given ε ∈ (0, 1), define g(ζ ) = exp[log(ε)(−q(ζ ))/(1 + M)]. Then g Ω ≤ 1 and g E ≤ ε, but |g(z)| = ε M/(M+1) . This contradicts (5.2). ⊓ ⊔ Together, Theorems 5.2 and 5.3 assert that analytic continuation with a boundedness condition in a nontrivial region is always well-posed but always has an infinite condition number. I am not aware if such a general assertion has been made before.

Analytic continuation in Chebfun
This project sprang from work with Chebfun, a software system for numerical computing with functions [9]. In its basic mode of operation, Chebfun works with smooth functions on an interval that without loss of generality we may take to be [−1, 1]. Chebfun represents each function to ≈ 16 digit precision by a polynomial in the form of a finite Chebyshev series, and early in the project, it was realized that the same series could be used for evaluation at complex points off the interval. Figure 6.1 illustrates this effect for the function f (x) = log(1 + x 2 ), which has branch points at x = ±i. By an adaptive process described in [1], the Chebfun command p = chebfun('log(1+x^2)') constructs a polynomial p that matches f on [−1, 1] with a maximal error of about 2 −51 ≈ 4.44 × 10 −16 ; the degree of p for this example is 38. Typing p(0) to evaluate the polynomial at x = 0, for example, returns the value 5.5 × 10 −17 , accurate to more than 16 digits. What is interesting is that typing p(i/2) also gives an accurate value: −0.2876820630, as compared with the true value log(0.75) ≈ −0.2876820768. How can we explain this?
In fact, Chebfun is carrying out the Chebyshev analogue of the algorithm of Section 3: it computes a finite sequence of Chebyshev series coefficients, then uses these coefficients to define an approximation p =f . Instead of working outward from the unit disk E to larger disks, it is working outward from the unit interval E to socalled Bernstein ellipses, whose algebra is defined by a transplantation of the results of Section 3 by the Joukowski map (z + z −1 )/2. For details, see [26], particularly the discussions of the "Chebfun ellipse" and the command plotregion. According to Theorem 3.1 as transplanted from disks to ellipses, we can expect the number of accurate digits to fall off smoothly, and Figure 6.1 shows that this is just what is observed. Analytic continuation in a region bounded by an ellipse is mentioned as Example 3 of [12], and a detailed analysis of algorithms in this geometry is presented in [7]. Results for ellipses analogous to those of [5] for disks can be found in [28].

Analytic continuation by a chain of disks
The classic idea for analytic continuation, going back to Weierstrass, involves a succession of Taylor expansions, each with its own disk of convergence. In principle, this procedure enables one to track a function along any path where it is analytic, and we recommend the beautifully illustrated discussion in Section 3.6 of Wegert's Visual Complex Functions [29]. For inexact function data, however, the method is far from promising. There is a small literature on numerical realizations, and a memorable contribution is a 1966 paper by Henrici in which the necessary transformations of series coefficients are formulated in terms of matrix multiplications; see [16] or [17, sec. 3.6].
To analyze this idea quantitatively, the simplest setting is a channel, essentially the same as the half-strip of Section 4. Specifically, consider the finite-length "stadium" G shown in Figure 7.1. Given L > 0, this is the strip of half-width 1 extending from x = 0 to x = L together with half-disks of radius 1 at each end. For some n > 0 and r ∈ (0, 1), we define h = L/n and for 0 ≤ k ≤ n. We assume n is large enough so that h < 1 − r. and then for all sufficiently small choices of ε, Proof We proceed by induction on k in (7.5). The case k = 1 follows from (7.3) by Theorem 3.1 (rescaled by a factor R = 1/r). Consider step k + 1, assuming (7.5) has been established for previous steps. Combining (7.4) and (7.5) gives By Theorem 3.1 (with the same rescaling as before), (7.6) implies (We avoid replacing log(1/r) by − log r since it can be confusing to have to remember that log r is negative.) We are done if we can show or by taking logarithms, By (7.2), if we divide both sides by the negative quantity log(ε k ), this becomes Now suppose for a moment that ε k is negligible. Then the log 2/ log(ε k ) term goes away and the condition we must verify reduces to A numerical search readily confirms that this holds with a strict inequality over the indicated region r ∈ (0, 1/2), h ∈ (0, 1/4) (the coefficient of the term −2h 2 was introduced to ensure this). Because of the assumption in the theorem statement that ε is sufficiently small, this establishes (7.8).

⊓ ⊔
The conclusion of Theorem 7.1 becomes memorable in the limit ε, h → 0. The parameterẽ of (7.2) is then minimized with the choice r = 1/e, for which it takes the valueẽ = e. We conclude that with this optimal choice of r, The number of accurate digits in chain-of-disks continuation along a channel of half-width 1 decays at the rate exp(−ex).
In a field as established as complex analysis, it is hard to be sure that anything is entirely new, but I am not aware of a previous estimation of loss of digits at the rate exp(−ex) for chain-of-disks continuation. The general observation of exponential loss of information is more than a century old. Henrici [16] writes (his italics) The early vectors. . . must be computed more accurately than the late ones and he gives credit for related work to Mittag-Leffler, Painlevé, Zeller, and Lewis [19].
It is interesting to note what our estimates suggest for what might be considered a very natural test problem for analytic continuation. Suppose we have a function like f (z) = √ z that is known to be known to be bounded and analytically continuable along any curve in the punctured disk 0 < |z| < 2. If we start near z = 1 with a certain accuracy ε and go around the origin and back to z = 1 again, how much accuracy will remain? We will not attempt to give a sharp solution to this problem, but it is the example that motivated our choice of a strip of length 2π for the numbers quoted in the introduction. Our estimates suggest that the number of accurate digits may be reduced by a factor as great as (π/4) exp(π 2 ) ≈ 15,000, or exp(2πe) ≈ 26,000,000 for the chain-of-disks method.

Conclusion
A compelling presentation of the practical side of analytic continuation can be found in the book to appear by Fornberg and Piret [11]. In the case of exactly known functions, although one could use the chain-of-disks idea in principle, Taylor series play little role in practice. A far more powerful approach is to find an analytical method to transform one formula defining a function (a formula being after all a finite object, unlike an infinite set of Taylor coefficients) into another formula with a new region of validity. For the most famous of all examples, the Dirichlet series for the Riemann zeta function converges only for Re z > 1, but other representations extend ζ (z) to the whole complex plane.
The present paper has concerned the case of inexactly known functions. Here, for continuation of functions from a disk to a larger disk, or from an interval to an ellipse, algorithms related to Taylor or Chebyshev series are effective, as has been discovered by various authors and we have illustrated by Figure 6.1 from Chebfun. A third equivalent context would be analytic continuation of a periodic function into a strip by Fourier series. The question is, what can one do to continue a function beyond the disk/ellipse/strip of convergence of its Taylor/Chebyshev/Fourier series? Our theorems show that if all one knows is analyticity and boundedness along certain channels, then accuracy may be lost at a precipitous exponential rate, better than the chain-of-disks but only by a constant. However, the assumption of analyticity just in a channel is more pessimistic than necessary in many applications. Functions arising in applications rarely have natural boundaries or other beautiful pathologies of analytic function theory, however generic such structures may be from a certain abstract point of view; they are far more likely to be analytic everywhere apart from certain poles and branch points. In practice, rational functions are the crucial tool for analytic continuation in such cases, and when they work, their convergence is typically exponential, just as we have found for series-based methods in a disk [3,10,26]. It would be an interesting challenge to develop theorems for meromorphic functions analogous to what we have established here in the analytic case, and a discussion with some of this flavor can be found in [21].