Printer Friendly


AMS subject classifications. 33C45, 65D15

1. Introduction. Let [mu] be a probability measure on R such that a family of [mathematical expression not reproducible]-orthonormal polynomials [mathematical expression not reproducible] can be defined.1 The non-decreasing function

[mathematical expression not reproducible]

is a probability distribution function on R since [mathematical expression not reproducible] has unit [mu]-integral over R. This paper is chiefly concerned with developing algorithms for drawing random samples of a random variable whose cumulative distribution function is [F.sub.n]. The high-level algorithmic idea is straightforward: develop robust algorithms for evaluating [F.sub.n], and subsequently use a standard root-finding approach to compute [F.sub.n.sup.-1]U) where U is a continuous uniform random variable on [0,1]. (This is colloquially called "inverse transform sampling".) The challenge that this paper addresses is in the computational evaluation of [F.sub.n] (x) for any n e [N.sub.0] and for relatively general [mu]. Borrowing terminology from [8], we call [F.sub.n] the order-n distribution induced by [mu].

In our algorithmic development, we focus on three classes of continuous measures [mu] from which induced distributions spring: (1) Jacobi distributions on [-1,1], (2) Freud (i.e., exponential) distributions on R, and (3) "Half-line" Freud distributions on [0, [infinity]). These measures encompass a relatively broad selection of continuous measures [mu] on R, including all measures associated with classical orthogonal polynomial families. Our algorithm applies directly to the measure classes identified in Table 1.1.

The utility of sampling from univariate induced distributions has recently come into light: The authors in various papers [11, 19, 3] note that additive mixtures of induced distributions are optimal sampling distributions for constructing multivariate polynomial approximations of functions using weighted discrete least-squares from independent and identically-distributed random samples. "Optimal" means that these distributions define a sampling strategy which

Freud [mathematical expression not reproducible] x[member of]U [rho ] >-1

[mathematical expression not reproducible] provides stability and accuracy guarantees with a sample complexity that is currently thought to be the best (smallest). This distribution also arises in related settings [12]. The ability to sample from an induced distribution, which this paper addresses, therefore has significant importance for multivariate applied approximation problems.

Induced distributions can also help provide insight for more theoretical problems. The weighted pluripotential equilibrium measure is a multivariate probability measure that describes asymptotic distributions of optimal sampling points [1, 19]. However, an explicit form for this measure is not known in general. The authors in [19] make conjectures about the Lebesgue density associated to equilibrium measure in one case when its explicit form is currently unknown. While these conjectures remain unproven, univariate induced distributions can be used to simulate samples from the equilibrium measure. Hence, sampling from induced distributions can be used to provide supporting evidence for the theoretical conjectures in [19].

The outline of this paper is as follows: in Section 2 we review many standard properties of general orthogonal polynomial systems that are exploited for computing induced distributions. Section 3 contains a detailed discussion of our novel approach for computing Fn(x) for three classes of measures; this section also utilizes potential theory results in order to approximate [F.sub.n.sup.-1](0.5). Section 4 uses the previous section's algorithms in order to formulate an algorithmic strategy for computing [F.sub.n.sup.-1], (u) [member of] [0,1]. Finally, Section 5 discusses the above-mentioned applications of multivariate polynomial approximation using discrete least-squares, and investigating conjectures for a weighted equilibrium measure.

Code that reproduces many of the plots in this paper is available for download [17]. The code contains routines for accomplishing almost all of the procedures in this paper including evaluation and inversion of induced distributions (for many of the distributions in Table 1.1), inverse transform sampling for multivariate sampling from additive mixtures of induced distributions, and fast versions of all codes that utilize approximate monotone spline interpolants for fast evaluation and inversion of distribution functions. The code also contains routines that reproduce Figure 1.1 (left, center), Figure 2.1 (right), Figure 3.1 (center, right), Figure 3.2 (left), Figure 3.3 (left), and Figure 5.1 (left).

1.1. A simple example. The main algorithmic novelties of this paper revolve around evaluation of [F.sub.n](x). In Figure 1.1 we show one example of the integrand [mathematical expression not reproducible](x)d[mu](x) and the associated [F.sub.n]. One suspects that packaged integration routines should be able to perform relatively well in order to compute integrals for such a problem. In our experience this is frequently true, but comes at a price of increased computational effort and time, and decreased robustness. The right-hand pane in Figure 1.1 shows timings for Matlab's built-in integral routine versus the algorithms developed in this paper. We see that the algorithms in this paper are much faster, usually resulting in around an order of magnitude savings when compared against integral. We caution that such tests compare an adaptive numerical integration method, Matlab's integral, against a tailored and optimized integration method, the algorithm of this paper. Therefore, an adaptive version of our algorithm is likely to be less efficient than reported here, and the non-adaptive version in this paper is naturally very specialized and therefore less useful for general integration problems.

Our experimentation (using Matlab's integral) also reveals the following advantages of using the specialized algorithm in this paper:

* The results from our algorithm appear to be more accurate compared to standard routines, and use significantly less computation. This conclusion is based on our testing with integral, and is true even if one modifies algorithm tolerances in integral.

* We have occasionally observed integral return non-monotonic evaluations for the distribution [F.sub.n]. This typically happens when most of the mass of the integrand is concentrated far from the boundary of the support of [mu]. Non-monotonic behavior causes problems in performing inverse transform sampling.

* When d[mu](x) has a singularity at boundaries of the support of [mu], then integral frequently complains about singularities and failure to achieve error tolerances. Similar results and conclusions hold when the algorithm of this paper is tested against Matlab's quadgk routine, another built-in integration function. We also note that the integrand [mathematical expression not reproducible](t)d[mu](t) can be highly oscillatory for large n, which suggests that asymptotics for this integrand along with techniques from oscillatory quadrature can be effective in computing these integrals. This paper does not investigate such approaches and instead focuses on exploiting properties of orthogonal polynomials. However, we believe that tools for oscillatory quadrature can be effective for this problem.

1.2. Historical discussion. The distribution function of the arcsine or "Chebyshev" measure is

[mathematical expression not reproducible]

It was shown in [20] that if [mu] belongs to the Nevai class of measures, then [F.sub.n](x) [right arrow] F(x) pointwise for all x e [-1,1]. Further refinements on this statement were made in [23] which generalized the class of measures for which convergence holds. Generating polynomials orthogonal with respect to the measure associated to [F.sub.n] is considered in [8] with a generalization given in [13].

The authors in [11] proposed sampling from an additive mixture of induced distributions using a Markov Chain Monte Carlo method for the purposes of computing polynomial approximations of functions via discrete least-squares; the work in [19] investigates sampling from the n-asymptotic limit of these additive mixtures. The authors in [3] leverage the additive mixture property to sample from this distribution using a monotone spline interpolant. At the very least this latter method requires multiple evaluations of [F.sub.n](x). To our knowledge there has been essentially no investigation into robust algorithms for the evaluation of [F.sub.n] for broad classes of measures, which is the subject of this paper.

2. Background.

2.1. Orthogonal polynomials. This section contains classical knowledge, most of which is available from any seminal reference on orthogonal polynomials [24, 4, 20, 6].

Let [mu] be a Lebesgue-Stieltjes probability measure on R, i.e., the distribution function

[mathematical expression not reproducible]

is non-decreasing and right-continuous on R, with F(-[infinity]) = 0 and F([infinity]) = 1. For any distribution F we use the notation [F.sup.c](x) := 1 - F(x) for its complementary function.

We assume that [mu] has an infinite number of points of increase, and has finite polynomial moments of all orders, i.e.,

[mathematical expression not reproducible]

Under these assumptions, a sequence of orthonormal polynomials [mathematical expression not reproducible] exists, with deg [p.sub.j] = j, satisfying

[mathematical expression not reproducible]

where [delta]k,j is the Kronecker delta function. We will write [p.sub.n](*) = [p.sub.n](*; [mu]) to denote explicit dependence of [p.sub.n] on [mu] when necessary. Such a family can be mechanically generated by iterative application of a three-term recurrence relation:

(2.1) [mathematical expression not reproducible]

where the recurrence coefficients an and bn are functions of the moments of [mu]. The initial conditions p-1 = 0 andp0 = 1 are used to seed the recurrence. With [p.sub.n] defined in this way, the (positive) leading coefficient of [p.sub.n] has value

[mathematical expression not reproducible]

The polynomial [p.sub.n] has n real-valued, distinct roots; when [mu] has a Lebesgue density with support given by a finite union of intervals, then the roots lie inside the convex hull of the support of [mu]. These roots [{[x.sub.k]}.sub.k=1.sup.n] are nodes for the Gaussian quadrature rule,

[mathematical expression not reproducible]

where the weights wk are unique and positive. These nodes and weights can be computed having knowledge only of the recurrence coefficients [a.sub.k] and [b.sub.k]; numerous modern algorithms accomplish this, with a historically significant procedure given in [10].

2.2. Induced orthogonal polynomials and measures. With [mathematical expression not reproducible] the orthonormal polynomial family with respect to [mu], the collection of polynomials orthogonal with respect to the weighted distribution [mathematical expression not reproducible] (x) with n fixed are called induced orthogonal polynomials. We adopt this terminology from [8].

Define the Lebesgue-Stieltjes measure [[mu].sub.n] and its associated distribution function [F.sub.n] by

(2.2) [mathematical expression not reproducible]

with [{pn{*)}n = {pn{*; [mu])}.sub.n] the orthonormal family for [mu]. Note that [mu]n [much less than] [mu], and Fn(co) = [[mu].sub.n] (R) = 1 so that [[mu].sub.n] is also a probability measure. The measure [[mu].sub.n] has its own three-term recurrence coefficients [a.sub.j,n] and [b.sub.j,n] for j = 0,..., that define a new set of [mathematical expression not reproducible]-orthonormal polynomials, which can be generated through the corresponding version of the mechanical procedure (2.1). One such procedure for generating these coefficients is given in [8].

We will call the measure [[mu].sub.n] the (order-n) induced measure for [mu], and [F.sub.n] the corresponding (order-n) induced distribution function.

Our main computational goal is, given u e [0,1], the evaluation of [mathematical expression not reproducible] (u) for various measures [mu]. The overall algorithm for accomplishing this is a root-finding method, e.g., bisection or Newton's method. Thus, the goal of finding [mathematical expression not reproducible] (u) also involves the evaluation of [F.sub.n](x), which is the focus of Section 3. A good root-finding algorithm also requires a reasonable initial guess for the solution. This initial guess is provided by the methodology in Section 4.1.

2.3. Measure modifications. Our algorithms rely on the ability to compute polynomial measure modifications. That is, given the three-term coefficients [a.sub.n] and [b.sub.n] for [mu], to compute the coefficients [a.sub.n] and [b.sub.n] for [mu] defined as

d[mu](x) = p(x)d[mu](x),

where p(x) is a polynomial, non-negative on the support of [mu]. This is a well-studied problem [9, 5, 6, 18]. In particular, one may reduce the problem to iterating over modifications by linear and quadratic polynomials. We describe in detail how to accomplish linear and quadratic modifications in the appendix, with the particular goal of structuring computations to avoid numerical under- and over-flow when n is large.

The following computational tasks described in the Appendix are used to accomplish measure modifications.

1. (Appendix A.1) Evaluation of [r.sub.n], the ratio of successive polynomials in the orthogonal sequence:

[mathematical expression not reproducible]

Above, we require x to lie outside the zero set of [p.sub.j-1].

2. (Appendix A.2) Evaluation of a normalized or weighted degree-n polynomial:

[mathematical expression not reproducible]

Note that [C.sub.n(x)/[r.sub.n](x) ~ 1 for large enough \x\.

3. (Appendix B) Polynomial measure modifications: given [mu] and its associated three-term recurrence coefficients [a.sub.n] and [b.sub.n], computation of the three-term recurrence coefficients associated with the measures [mu] and [mu], defined as d[mu](x) = [+or -] (x - [y.sub.0]) d[mu](x), [y.sub.0] [member of] supp[mu], d[mu]{x) = (x - z0) d[mu](x), [z.sub.0] [member of] R.

The + sign in [mu] is chosen so that [mu] is a positive measure.

3. Evaluation of [F.sub.n]. This section develops computational algorithms for the evaluation of the induced distribution [F.sub.n] defined in (2.2). These algorithms depend fairly heavily on the form of a Lebesgue density d[mu](x) (i.e., a positive weight function) for the measure [mu] on R, but the ideas can be generalized to various measures. We consider the classes of weights enumerated in Table 1.1:

* (Jacobi weights) d[mu]J(x) = (1 - x)[alpha](1 + x)[beta] for x [member of] [-1,1] with parameters [alpha],[beta] > -1.

* (Freud weights) d[mu]F(x) = |x|p exp(- |x|[alpha]) for x [member of] R with parameters [alpha] > 0, [rho ] >-1.

* (half-line Freud weights) d[mu]HF(x) = [x.sup.p] exp (-[x.sup.a]) for x e [0, [infinity]) with parameters [alpha] > 0, [rho ] > -1.

The induced distribution Fn for Freud weights can actually be written explicitly in terms of the corresponding induced distribution for half-line Freud weights, so most of the algorithm development concentrates on the Jacobi and half-line Freud cases. The strategy for these two latter cases is essentially the same: with [mu] in one of the classes above and n fixed, we divide the computation into one of two approximations, depending on the value of x. Each approximation is accurate for its corresponding values of x. Formally, the algorithm is

(31) [mathematical expression not reproducible]

where [F.sub.n](x) and [mathematical expression not reproducible] represent computational approximations to [F.sub.n](x) and [mathematical expression not reproducible], respectively, and are the outputs from the algorithms that we will develop. A pictorial description of this is given in Figure 2.1 (left). We choose to seek a value of [x.sub.0] such that [mathematical expression not reproducible]. This is a heuristic choice that we observe works well in practice across all the measures we consider. It is likely that other choices for [x.sub.0] can also produce effective algorithms, but we have found our choice to be simple and amenable to direct computation. However, we cannot know this value [mathematical expression not reproducible] (0.5) a priori, so we use potential-theoretic arguments to compute a value [x.sub.0] = [x.sub.0] p [mu], n) approximating the median of [F.sub.n],

(3.2) [mathematical expression not reproducible]

For our choice of [x.sub.0], we can provide no estimates for this approximate equality, but empirical evidence in Figure 2.1 (right) shows that our choices are very close to the real median, uniformly in n. The coming sections concentrate on, for each class of [mu] mentioned above, specifying [x.sub.0] and detailing algorithms for [F.sub.n] and [mathematical expression not reproducible].

3.1. Jacobi weights. We consider computing induced distribution functions for Jacobi

measures [mathematical expression not reproducible] as defined in Table 1.1. When circumstances are clear, we will write [mathematical expression not reproducible] to avoid notational clutter. We seek the distribution function of the induced measure,

(3.3) [mathematical expression not reproducible]

To compute [F.sub.n], we specify an approximate median [x.sub.0] = [x.sub.0] p [mu], n) satisfying (3.2), and construct algorithmic procedures to evaluate [E.sub.n](x) [much less than] [F.sub.n](x) (for x [??] [x.sub.0]) and [mathematical expression not reproducible] (for x > [x.sub.0]). Having specified these, we use (3.1) to compute our approximation to [F.sub.n].

3.1.1. Computation of [x.sub.0]pn). With [alpha], [beta] > - 1 and n P IN fixed, consider the measure [mathematical expression not reproducible]. Note that this measure is still a Jacobi measure since [mathematical expression not reproducible] > - 1 and [mathematical expression not reproducible] > -1. Since

[mathematical expression not reproducible]

the integrand in (3.3) satisfies

[mathematical expression not reproducible]

The quantity under the square brackets on the right-hand side is, in the language of potential theory, a weighted polynomial of degree n. One result in potential theory characterizes the "essential" support of this weighted polynomial; in particular, the weighted polynomial decays quickly outside this support. The essential support is an interval, and we take the median [x.sub.0] of the induced measure [[mu].sub.n] to be the centroid of this interval.

The essential support of the weighted polynomial above is demarcated by the Mhaskar-Rakhmanov-Saff numbers for the asymmetric weight [mathematical expression not reproducible] on [-1,1]. Thesenumbers for this weight are computed explicitly in [22, pp 206-207]. When [alpha],[beta] [??] 0, the support interval is [[theta] - [DELTA], [theta] + [DELTA]], with

[mathematical expression not reproducible]

Since this is the interval where most of the "mass" of the integral in (3.3) lies, we set [x.sub.0] to be the centroid [theta] of this interval:

(3.4) [mathematical expression not reproducible]

Note that the definitions of [theta] and [DELTA] require [alpha] and [beta] to be non-negative. Without mathematical justification, we extend usage of the formula (3.4) to valid negative values of [alpha], [beta] as well. Figure 2.1 compares [x.sub.0] and F~1 (0.5) for certain choices of [alpha] and [beta].

3.1.2. Computation of [F.sub.n](x). First assume that x [??] [x.sub.0], with [x.sub.0] defined in (3.4), and define A [member of] [N.sub.0] as

A:=||[alpha]| [right arrow] [alpha] - A [member of] (-1,1),

where [*] is the floor function. We transform the integral (3.3) over [- 1, x] onto the standard interval [- 1,1] via the substitution [mathematical expression not reproducible]:

[mathematical expression not reproducible]

where [U.sub.2n+A](u) is a degree-(2n + A) polynomial given by (3.5) [mathematical expression not reproducible]

[mathematical expression not reproducible]

where [mathematical expression not reproducible] are the n zeros of pn(*). Since we explicitly know the polynomial roots of [U.sub.2n+A], we can absorb the factor marked (a) into the measure [d[mu].sup. (0,[beta])(u) via A linear modifications, and we can likewise absorb the factor (b) via n quadratic modifications. (See Appendix B.) Thus, define

(3.6) d[[mu].sub.n](u) = [U.sub.2n+A](u)[d[mu].sub.(0,[beta])(u),

which is a modified measure whose recurrence coefficients can be computed via successive application of the linear and quadratic modification methods in Appendix B. Thus,

(3.7) [mathematical expression not reproducible]

[mathematical expression not reproducible]

The integrand above has a root ([alpha] > 0) or singularity ([alpha] < 0) at [mathematical expression not reproducible] this root is far outside the interval [-1,1] unless [beta] is very large and both n and [alpha] are small. The integrand is therefore a positive, monotonic, smooth function on [- 1,1], taking values between 1 - x and 2; we use an order-M[[mu]n]-Gaussian quadrature to efficiently evaluate it. With [mathematical expression not reproducible] denoting the nodes and weights, respectively, of this quadrature rule, we compute

(3.8) [mathematical expression not reproducible]

(3.9) [mathematical expression not reproducible]

The entire procedure is summarized in Algorithm 1.

In order to compute [mathematical expression not reproducible] for x > [x.sub.0], we use symmetry. Since x [right arrow]x interchanges the parameters [alpha] and [beta], then

[mathematical expression not reproducible]

Note that if x > [x.sub.0] ([[mu].sup.([alpha],[beta])] n), then - x < [x.sub.0] ([[mu].sup.([beta],[alpha])] n). Thus, [mathematical expression not reproducible] can be computed via the same algorithm for [F.sup.n], but with different values for [alpha], [beta], and x. This is also shown in Algorithm 1.

Algorithm 1: Computation of [mathematical expression not reproducible], approximating [F.sub.n](x) for [mu] corresponding to aJacobi polynomial measure.

input : [alpha], [beta] > - 1: Jacobi polynomial parameters.

input :n [member of] [N.sub.0] and x [member of] [-1,1]: Order of induced polynomial and measure [[mu].sub.n] and value x.

input : M [member of] N: Quadrature order for approximate computation of [F.sub.n](x).

output :[F.sub.n](x;[[mu].sup.([alpha],[beta])]).

1 If x > [x.sub.0] {[mu]{[alpha],[beta]),n), return 1 - [F.sub.n](-x; [mu]([beta],[alpha])).

2 Compute n zeros, [mathematical expression not reproducible] of pn = pn (*; [[mu].sup.([alpha],[beta])]), and leading coefficient [[gamma].sub.n] of [p.sub.n].

3 Compute recurrence coefficients [a.sub.j] and [b.sub.j] associated to [[mu].sub.(0,[beta])] for 0 [??] j [??] M + A + 2n.

4 for j = 1,... ,n do

5 Quadratic measure modification (B.1b): update an, bn for

n = 0,..., M + A + 2(n - j) with modification factor

[mathematical expression not reproducible]

6 end

7 for k = 1,..., Ado

8 Use linear modification (B.1a) to update [a.sub.n], [b.sub.n] for n = 0,..., M + (A - k) with modification factor [mathematical expression not reproducible].

9 end

10 [mathematical expression not reproducible] (Cf. the leading coefficient of the polynomial [U.sub.2n+A] in (3.5).)

11 Compute M-point Gauss quadrature [mathematical expression not reproducible] associated with measure [[mu].sub.n] via [mathematical expression not reproducible].

12 Compute the integral I in (3.8), and return [mathematical expression not reproducible] given by (3.9).

THEOREM 3.1. With [[mu].sup.([alpha],[beta]) n [member of] N, and x [member of] [-1,1] all given, assume that x [??] [x.sub.0] with [x.sub.0] as in (3.4). Then the output [F.sub.n]{x) from Algorithm 1 using an M-point quadrature rule (3.10) [mathematical expression not reproducible]

where [b.sub.j] ([[mu].sub.n]) are the [b.sub.j] three-term recurrence coefficients associated to the x-dependent measure [[mu].sub.n] defined in (3.6). The constant C is

[mathematical expression not reproducible]

Note that C([alpha], [beta], n, M) on the right-hand side of (3.10) is explicitly computable once [alpha], [beta], M, and n are fixed, and only the product involving bj ([[mu].sub.n]) depends on x. Furthermore, the [b.sub.j] ([[mu].sub.n]) factors are explicitly computed in Algorithm 1, so that a rigorous error estimate for the algorithm can be computed before its termination.

Since [mathematical expression not reproducible] < 1/2, then the estimate (3.10) also hints at exponential convergence of the quadrature rule strategy, assuming that the factors [b.sub.j] ([[mu].sub.n]) can be bounded or controlled. We cannot provide these bounds, although we do know the asymptotic behavior [b.sub.j] ([[mu].sub.n]) [right arrow] 1/4 as j [right arrow] [infinity] [20, Remark 3.1.10]. Also, since [x.sub.0] e [-1,1] for any n > 0 then, uniformly in n > 0, C([alpha], [beta], n, M) [??] C'([alpha], [beta])4-M. Thus, for a fixed x we expect that the estimate in Theorem 3.1 behaves like

[mathematical expression not reproducible]

showing exponential convergence with respect to M. However, we cannot prove this latter statement.

Proof of Theorem 3.1. The result is a relatively straightforward application of known error estimates for Gaussian quadrature with respect to non-classical weights. We use the notation of Algorithm 1: [u.sub.m] and [w.sub.m] denote the M-point [[mu].sub.n]-Gaussian quadrature nodes and weights, respectively. We start with the Corollary to Theorem 1.48 in [6], stating that if f(*) is infinitely differentiable on [-1,1], then

[mathematical expression not reproducible]

for some [tau] [member of] (-1,1). Noting that since [u.sub.m] are the zeros of the degree-M [[mu].sub.n]-orthogonal polynomial, then

[mathematical expression not reproducible]

From (3.7), the integral that we would like to approximate has the integrand

[f(u) = (2 - 1/2(u + 1)(x + 1)).sup.[alpha]-A] . Then for any [tau] e (-1,1) and any x [??] [x.sub.0],

[mathematical expression not reproducible]


[mathematical expression not reproducible]

Since [alpha] - A [member of] (-1, 1), then

[mathematical expression not reproducible]

The result follows by direct computation of [mathematical expression not reproducible]

We verify spectral convergence of the scheme empirically in Figure 3.1, which also illustrates that for the test cases shown, one could choose M to bound errors uniformly in x.

The figure also shows that qualitative error behavior is uniform even for extremely large values of n and/or a or [beta]. Based on these results, taking M = 10 appears sufficient uniformly over all (a, [beta], n). (2) Finally, extending the estimate in (3.10) to the case x > [x.sub.0] can be accomplished by permuting a and [beta] as is done for that case in Algorithm 1.

3.2. Half-line Freud weights. In this section we consider the half-line Freud measure [mathematical expression not reproducible] defined in Table 1.1. These algorithms require the recurrence coefficients for [mathematical expression not reproducible] these coefficients are in general not easy to compute when a [not equal to] 1. We show in Appendix C that these recurrence coefficients can be determined from the recurrence coefficients for Freud weights, but recurrence coefficients for Freud weights are themselves relatively difficult to tabulate [15, 25]. In our computations we use the methodology of [7] to compute Freud weight recurrence coefficients (and hence half-line Freud coefficients). We note that the methodology of [7] is computationally onerous: For a fixed a and p, it required a day-long computation to obtain 500 recurrence coefficients.

Again in this subsection we write [mathematical expression not reproducible] similarly for [u.sub.n], F, [F.sub.n], etc. We restore these super- and subscripts when ambiguity arises without them.

We accomplish computation of Fn for this measure with largely the same procedure as for Jacobi measures. Like in the Jacobi case, the details of the procedure we use differ depending on whether x is closer to the left-hand end of supp [mu] (which is x = 0 here), or to the right-hand end of the supp [mu] (which is x = [infinity]). We determine this delineation again by means of potential theory.

3.2.1. Computation of [X.sub.0]. As with the Jacobi case, we take [x.sub.0] to be the midpoint of the "essential" support for [mathematical expression not reproducible], the latter of which is approximately the support of the weighted equilibrium measure associated to [mathematical expression not reproducible]. However, directly computing the support of this equilibrium measure is difficult. Thus, we resort to a more ad hoc approach.

To derive our approach, we first compute the exact support of special cases of our Half-line Freud measures.

* The support of the weighted equilibrium measure associated to the measure

(3.11) d[mu](x) = [x.sup.s] exp(-[lambda]x), x [member of] [0, [infinity]), s [??] 0, [lambda] > 0,

is the interval [[theta] - [DELTA], [theta] + [DELTA]], with these values given by [22]

[mathematical expression not reproducible]

This interval is the essential support for any function of the form pn(x) (d[mu](x)) where [p.sub.n] is a degree-n polynomial. * The second special case is for arbitrary [alpha], but [rho ] = 0. The support of the weighted equilibrium measure for [mathematical expression not reproducible] in this case is the interval [0, [k.sub.n]([alpha])], where

[mathematical expression not reproducible]

are the Mhaskar-Rakhmanov-Saff numbers for [mathematical expression not reproducible] (see [16]).

We now derive our approximation for the case of general [alpha], [rho ] . To approximate where [mathematical expression not reproducible] is supported, consider

[mathematical expression not reproducible]

where we have introduced [q.sub.n/a], which is a "polynomial" of "degree" n/[alpha] . (3) Note that in the variable u, the weight function under square brackets in the last expression is of the form (3.11). Concepts in potential theory extend to generalized notions of polynomial degree, and so we may apply our formulas for [theta] and [DELTA] with s = p/2n and [lambda] = [alpha]/2n. These formulas imply that the "essential" support for the u variable is

[theta] - [DELTA] [??] u = [x.sup.[alpha] [??] [theta] + [DELTA].

Therefore, to obtain appropriate limits on the variable x, we raise the endpoints [theta] + [DELTA] to the 1/[alpha] power. However, we now require a correction factor. To see why, we compute the right-hand side of our computed support interval when [rho ] = 0

(3.12) [mathematical expression not reproducible]

and compare this with the exact value kn ([alpha]) computed above. We note that while [k.sup.n] ([alpha]) ~ [n.sup.1la] matches the n-behavior of (3.12), the constant is wrong. We thus multiply the endpoints [([theta] + [DELTA]).sup.1/a] by the appropriate constant to match the [rho ] = 0 behavior of kn([alpha]). The net result then, for arbitrary [alpha], [rho ] , is the approximation

(3.13a) [mathematical expression not reproducible]

(3.13b) [mathematical expression not reproducible]

Figure 3.2 compares the intervals demarcated by a_ and a+ versus [mathematical expression not reproducible] ([0.01,0.99]), the latter of which contains "most" of the support of [F.sub.n].

3.2.2. Computation of [F.sub.n](x). First assume that x [??] [x.sub.Q]. Then

[mathematical expression not reproducible]t

We recognize a portion of the integrand as a Jacobi measure, and use successive measure modifications to define [[mu].sub.n]:

[mathematical expression not reproducible]

(3.14a) [mathematical expression not reproducible]

(3.14b) [mathematical expression not reproducible]

The recurrence coefficients of [[mu].sub.n] can be computed via polynomial measure modifications on the roots [mathematical expression not reproducible], where [x.sub.j,n] the roots of [p.sub.n] (*). A detailed algorithm is given in Algorithm 2.

Algorithm 2: Computation of [F.sub.n] (x) for [mathematical expression not reproducible] corresponding to a half-line Freud weight.

input :[alpha] > 1/2, [rho ] > - 1: half-line Freud weight parameters.

input : n [member of] [N.sub.0] and x [??] 0: Order of induced measure [[mu].sub.n] and value x.

input : M [member of] N: Quadrature order for approximate computation of [F.sub.n](x). output :F n(x).

1 Compute n zeros, [mathematical expression not reproducible] of pn = [p.sub.n] (*; [[mu].sup.(a,[rho ] )]), and leading coefficient [[gamma]n] of pn.

2 Compute recurrence coefficients [a.sub.j] and [b.sub.j] associated to [mathematical expression not reproducible] for 0 [??] j [??] M + 2n.

3 for j = 1,..., n do

4 Quadratic measure modification (B.1b): update [a.sub.n] and [b.sub.n] for n = 0,..., M + 2(n - j) with modification factor [mathematical expression not reproducible].

5 Scale [b.sub.0] [left arrow] [b.sub.0 exp] (1/n log [[gamma].sub.n/2).

6 end

7 Compute M-point Gauss quadrature [mathematical expression not reproducible] associated with measure [[mu].sub.n] via [mathematical expression not reproducible]

8 Compute the integral I in (3.14a), and return [mathematical expression not reproducible] given by (3.14b).

Now assume that x [??] [x.sub.0]. We compute [mathematical expression not reproducible] directly in a similar fashion as we did for [F.sub.n]. We have

[mathematical expression not reproducible]

We again use this to define a new measure [[mu].sub.n] and an associated M-point Gauss quadrature (um, wm)m=1. The recurrence coefficients for [[mu].sub.n] are computable via polynomial measure modifications. This results in the approximation

[mathematical expression not reproducible]

(3.15) [mathematical expression not reproducible]

(3.16) [mathematical expression not reproducible]

A more detailed algorithm is given in Algorithm 3. Of course, once [mathematical expression not reproducible] is computed we may compute [mathematical expression not reproducible].

Algorithm 3: Computation of [mathematical expression not reproducible] for [mathematical expression not reproducible] corresponding to a half-line Freud weight.

input :[alpha] > 1/2 [rho ] > - 1: generalized Freud weight parameters.

input :n [member of] [N.sub.0] and x [??] 0: Order of induced polynomial and measure [[mu].sub.n] and value x.

input : M e IN: Quadrature order for approximate computation of Fc (x). output :Fc{x).

1 Compute n zeros, [mathematical expression not reproducible] of [mathematical expression not reproducible], and leading coefficient [[gamma].sub.n] ofpn.

2 Compute recurrence coefficients aj and bj associated to [mathematical expression not reproducible] for 0 [??] j [??] M + 2n.

3 for j = 1,..., n do

4 Quadratic measure modification (B.1b): update [a.sub.n], [b.sub.n] for n = 0,..., M + A + 2(n - j) with modification factor (u - ([x.sub.j,n] - x)) (2).

5 Scale [b.sub.0] <- [b.sub.0] exp (1/n log [mathematical expression not reproducible]).

6 end

7 Compute M-point Gauss quadrature [mathematical expression not reproducible] associated with measure [[mu].sub.n] via [mathematical expression not reproducible].

8 Compute the integral I in (3.15), and return [mathematical expression not reproducible] given by (3.16).

Differences between [Fn.sub.n] and computational approximations [F.sub.n] are shown in Figure 3.3. We see that we require a much larger value of M in order to achieve accurate approximations compared to the Jacobi case. We believe this to be the case due to the function exp(u[alpha] + x[alpha] -(u + x)[alpha]) appearing in the integral I(x). Note that for [alpha] = 1 this function becomes unity and so does not adversely affect the integral; this results in the much more favorable error plot on the left in Figure 3.3.

The different behavior for [alpha] = 1 leads us to make customized choices in this case: we choose M = 25 for all values of n and [rho ] , and we take x0(n) = 50. Whereas for [alpha] [not equal to] 1 our tests suggest that M = n + 10 is sufficient to achieve good accuracy, and we take x0 as the average of a+ as given in (3.13).

3.3. Freud weights. Finally, consider the Freud measure pF, defined in Table 1.1. An especially important case occurs for [alpha] = 2, [rho ] = 0, corresponding to the classical Hermite polynomials. Note that the recurrence coefficients for general values of a and [rho ] are not known explicitly, but their asymptotic behavior has been established [14].

It is well-known that Freud weights are essentially half-line Freud weights in disguise under a quadratic map. It is not surprising then that we may express primitives of induced polynomial measures for Freud weights in terms of the associated half-line Freud primitives.

THEOREM 3.2. Let the parameters (a, p) define a Freud weight and associated measure [mathematical expression not reproducible]. Then for x [??] 0,

(3.17) [mathematical expression not reproducible]

For x [??] 0, we have (3.18) [mathematical expression not reproducible]

Note that with expressions (3.17) and (3.18), an algorithm for computing [F.sub.n] (*; [[mu]F]) is straightforward to devise utilizing Algorithm 3 for [mathematical expression not reproducible] (*; [mu]HF).

The result (3.18) follows easily from the fact that the integrand in (2.2) defining [F.sub.n] is an even function. To prove the main portion of the theorem, expression (3.17), we require the following result relating Freud orthonormal polynomials to half-line Freud orthonormal polynomials.

LEMMA 3.3. Let p > - 1 and a > 1 be parameters that define a Freud measure [mathematical expression not reproducible] with associated orthonormal polynomial family [mathematical expression not reproducible]. Define two sets of half-line Freud parameters (a*, p*) and (a**, [rho]**) and the corresponding half-line Freud measures and polynomials:

(3.19a) [mathematical expression not reproducible]

(3.19b) [mathematical expression not reproducible]

Also define the constant

(3.20) [mathematical expression not reproducible]

Then, for all n [??] 0,

[mathematical expression not reproducible]

Proof. The following equalities may be verified via direct computation using the definitions (3.19) and (3.20) along with the expressions in Table 1.1:

(3.21) [mathematical expression not reproducible]

The proof of this lemma relies on the change of measure t [right arrow] [x.sup.2]. We have

[mathematical expression not reproducible]

This relation shows that the family [mathematical expression not reproducible] are polynomials of degree 2n that are orthonormal under a Freud weight with parameters [alpha] = 2[alpha]* and [rho ] = 2[rho ] * + 1. Using nearly the same arguments, but with the family {p**,n}n, yields the relation

[mathematical expression not reproducible]

This relation shows that the family [mathematical expression not reproducible] are polynomials of degree 2n + 1 that are orthonormal under the same Freud having parameters [alpha] = 2[alpha]** and [rho ] = 2[rho ] ** - 1. We also have that xp**,n ([x.sub.2]) is orthogonal to p*,m ([x.sub.2]) under a(ny) Freud weight for any n, m because of even-odd symmetry. Thus, define

Then [mathematical expression not reproducible] is a family of degree-n polynomials (with positive leading coefficient) orthonormal under a ([alpha], [rho ] ) Freud weight. Therefore, [P.sub.n] = [p.sub.n].

We can now give the proof of Theorem 3.2.

Proof of Theorem 3.2. Assume x [??] 0. Then

[mathematical expression not reproducible]

where ([alpha]*, [rho ] *) are as defined in (3.19a). By Lemma 3.3,

[mathematical expression not reproducible]

where ([alpha]**, [rho ] **) are defined in (3.19b). Then if n is even, we have

[mathematical expression not reproducible]

where we recall the equalities (3.21) related cF to cF*. Similarly, if n is odd we have

[mathematical expression not reproducible]

The combination of these results proves (3.17).

4. Inverting induced distributions. We have discussed at length in previous sections algorithms for computing [F.sub.n] px) defined in (2.2) for various Lebesgue-continuous measures [mu]. The central application of these algorithms we investigate in this paper is actually in the evaluation of [mathematical expression not reproducible] for u [member of] 0,1]. We accomplish this by solving for x in the equation

(4.1) [E.sub.n](x)-u = 0, u [member of] [0,1],

using a root-finding method. Our first step involves providing an initial guess for x.

4.1. Computing an initial interval. We use s+ to denote the (possibly infinite) end-points of the support of [mu] :

- [infinity] [??] s_ := inf (supp [mu]), s+ := sup (supp [mu]) [??] [INFINITY].

Now let u [member of] [0,1]. Our first step in finding [mathematical expression not reproducible] (u) is to compute two values x- and x+, such that

(4.2) [mathematical expression not reproducible]

Our procedure for identifying an initial interval containing [mathematical expression not reproducible] leverages the Markov-Stieltjes inequalities for orthogonal polynomials. These inequalities state that empirical probability distributions of Gauss quadrature rules generated from a measure bound the distribution function for this measure. Precisely:

LEMMA 4.1 (Markov-Stieltjes Inequalities, [24]). Let [mu] be a probability measure on R with an associated orthogonal polynomial family. For any N [member of] N, let [mathematical expression not reproducible] denote the N [mu]-Gaussian quadrature nodes and weights, respectively. Then:

[mathematical expression not reproducible]

Let [mathematical expression not reproducible] denote the sequence of polynomial orthonormal under the induced measure [[mu].sub.n]. Given N e IN, let [mathematical expression not reproducible] denote the N -point [[mu].sub.n]-Gaussian quadrature rule, i.e., [z.sub.k],N are the N ordered zeros of pN,n(x), with [v.sub.k],N the associated weights. Since [[mu].sub.n] is a probability measure, then [mathematical expression not reproducible], N = 1. As such, given u [member of] [0,1] we can always find some m [member of] {1,..., N}, such that

(4.3) [mathematical expression not reproducible]

Then, defining [z.sub.0],N = s- and [z.sub.N+1], N = s+ for all N and n, we have

[mathematical expression not reproducible]

Since [F.sub.n] is non-decreasing, this is equivalently,

[mathematical expression not reproducible].

Thus, if we find an m such that (4.3) holds, then (4.2) holds with

(4.4) [mathematical expression not reproducible].

When supp [mu] is bounded, the N -asymptotic density of orthogonal polynomial zeros on supp [mu] guarantees that we can find a bounding interval with endpoints x+ of arbitrarily small width by taking N sufficiently large. The difficulty is that we therefore require the zeros zk,N and the quadrature weights [v.sub.k],N of the induced measure, which in turn require knowledge of the three-term recurrence coefficients associated to [[mu].sub.n]. These can be easily computed from the coefficients associated to [mu]; since

(4.5) [mathematical expression not reproducible]

then we may again iteratively utilize the quadratic modification algorithm given by (B.1b) to compute these recurrence coefficients, which are iteratively quadratic modifications of the [mu]-coefficients. (Note that this is precisely the procedure proposed in [8] for computing these coefficients.)

4.2. Bisection. For simplicity, the root-finding method we employ to solve (4.1) is the bisection approach. More sophisticated methods may be used, with the caveat that the derivative of the function, [mathematical expression not reproducible], vanishes wherever [p.sub.n] has a root. We have found that a naive application of Newton's method for root-finding often runs into trouble, even with a very accurate initial guess.

The bisection method for root-finding applied to (4.1) starts with an initial guess for an interval [x-, x+] containing the root x, and iteratively updates this interval via

[mathematical expression not reproducible]

After a sufficient number of iterations so that x+ - x+ and/or |F(x-) - F(x+)| is smaller than a tunable tolerance parameter, one can confidently claim to have found the root x to within this tolerance. A good initial guess for x[+ or -] lessens the number of evaluations of [F.sub.n] in a bisection approach and thus accelerates overall evaluation of [mathematical expression not reproducible].

4.3. Algorithm summary. The overall algorithm for solving (4.1) is to (i) compute the recurrence coefficients associated with [[mu].sub.n] in (4.5) via quadratic measure modifications, (ii) compute order-N [[mu].sub.n]-Gaussian quadrature nodes and weights [z.sub.j,N] and [v.sub.j,N], respectively, (iii) identify m such that (4.3) holds so that x+ may be computed in (4.4), and (iv) iteratively apply the bisection algorithm with the initial interval defined by x+ using the evaluation procedures for [F.sub.n] outlined in Section 3.

5. Applications. This section discusses two applications of sampling from univariate induced measures. Both these applications consider multivariate scenarios, and are based on the fact that many "interesting" multivariate sampling measures are additive mixtures of tensorized univariate induced measures. Our first task is to introduce notation for tensorized orthogonal polynomials.

We will write d-variate polynomials using multi-index notation: [lambda] P [mathematical expression not reproducible] denotes a multi-index with components [lambda] = ([lambda]i,..., [lambda]d) and magnitude [mathematical expression not reproducible]. A point x [member of] [R.sup.d] has components x = ([x.sub.1],..., [x.sub.d]), and [mathematical expression not reproducible]. A collection of multi-indices will be denoted [mathematical expression not reproducible]; we will assume that N = |[LAMBDA]| is finite.

Let [mu] be a tensorial measure on Rd such that each of its d marginal univariate measures [[mu].sup.(j), j = 1,..., d admits a [[mu].sup.(j)-orthonormal polynomial family [mathematical expression not reproducible] on R, satisfying

[mathematical expression not reproducible]

A tensorial [mu] allows us to explicitly construct an orthonormal polynomial family for [mu] from univariate polynomials,

[mathematical expression not reproducible]

These polynomials are an [mathematical expression not reproducible]-orthonormal basis for the subspace PA, defined as

PA = span {p[lambda] | [lambda] [member of] [LAMBDA]}.

Under the additional assumption that the index set [LAMBDA] is downward-closed, then PA = span {x[lambda] | [lambda] [member of] [LAMBDA]}.

We extend our definition of induced polynomials to this tensorial multivariate situation. For any [lambda] P [LAMBDA], the order-[lambda] induced measure [mu][lambda] is defined as

[mathematical expression not reproducible]

where [mathematical expression not reproducible] is the (univariate) order-[lambda]j induced measure for [[mu].sup.(j)] according to the definition (2.2). Thus, [[mu].sub.[lambda]] is also a tensorial measure.

5.1. Optimal polynomial discrete least-squares. The goal of this section is to describe a procedure utilizing the algorithms above for performing discrete least-squares recovery in a polynomial subspace using the optimal (fewest) number of samples. The procedure we discuss was proposed in [3] and is based on the foundational matrix concentration estimates for least-squares derived in [2].

Let f : [R.sup.d] [right arrow] R be a d-variate function. Given (i) a tensorial probability measure [mu] admitting an orthonormal polynomial family, and (ii) a dimension-N polynomial subspace PA, we are interested in approximating the [mathematical expression not reproducible] -orthogonal projection of f onto PA. This projection is given explicitly by

[mathematical expression not reproducible]

One way to approximate the integral defining the coefficients [mathematical expression not reproducible] is via a Monte Carlo least-squares procedure using M collocation samples of the function f(x). Let [mathematical expression not reproducible] denote a collection of M independent and identically distributed random variables on [M.sup.d], where we leave the distribution of [X.sub.m] unspecified for the moment. A weighted discrete least-squares recovery procedure approximates [mathematical expression not reproducible] with c[lambda], computed as

(5.1) [mathematical expression not reproducible]

where [w.sub.m] are positive weights. One supposes that if the distribution of [X.sub.m] and the weights [w.sub.m] are chosen intelligently, then it is possible to recover the N coefficients [c.sub.[lambda]] with a relatively small number of samples M; ideally, M should be close to N. The analysis in [2] codifies conditions on a required sample count M so that the minimization procedure above is stable, and so that the recovered coefficients [c.sub.[lambda]] are "close" to [mathematical expression not reproducible]; these conditions depend on the distribution of [X.sub.m], on [w.sub.m], on [mu], and on PA. Since [mu] and [P.sub.A] are specified, the goal here is identification of an appropriate distribution for [X.sub.m] and weight [w.sub.m].

Using ideas proposed in [19, 11] the results in [3] show that, in the context of the analysis in [2], the optimal choice of probability measure [mu]X for sampling [X.sub.m] and weights [w.sub.m] that achieves a minimal sample count M are

(5.2) [mathematical expression not reproducible]

The precise quantification of the sample count and error estimates can be formulated using an algebraic characterization of (5.1). Define matrices V [member of] [R.sup.MxN] and W [member of] [K.sup.MxM], and vectors c e [R.sup.N] and f [member of] [R.sup.M] as follows:

[mathematical expression not reproducible]

where [lambda](1),..., [lambda](N) represents any enumeration of elements in [LAMBDA]. We use || * || on matrices to denote the induced [l.sup.2] norm. The algebraic version of (5.1) is then to compute c that minimizes the the least-squares residual of [mathematical expression not reproducible]. The following result holds.

THEOREM 5.1 ([2, 3]). Let 0 < [delta] < 1 and r > 0 be given, and define [c.sub.[delta]] := (1 + [delta]) log(1 + [delta]) - [delta] e (0,1). Draw M iid samples [mathematical expression not reproducible] from [mu]X, and let the coefficients c[lambda] be those recovered from (5.1). If

(5.3) [mathematical expression not reproducible]


[mathematical expression not reproducible]

The free parameter r is a tunable oversampling rate; [delta] represents the guaranteeable proximity of [V.sup.T]WV to I. We emphasize that by choosing [[mu].sub.X] = [[mu].sub.[LAMBDA]] with the weights defined as in W, then the size of M as required by (5.3) depends only on the the cardinality N of [LAMBDA], and not on its shape. Furthermore, the criterion M/ log M > N is optimal up to the logarithmic factor. Also, the statements above hold uniformly over all multivariate [mu].

Note that the optimal sampling measure [[mu].sub.X] is an additive mixture of induced measures and can be easily sampled, assuming that [[mu].sub.[lambda]] can be sampled. Sampling from [[mu].sub.X] defined above is fairly straightforward given the algorithms in this paper: (1) given [LAMBDA] choose an element [lambda] randomly using the uniform probability law, (2) generate d independent, uniform, continuous random variables Uj, j = 1,...,d, each on the interval [0,1], (3) compute X [member of] [R.sup.d] as

[mathematical expression not reproducible]

Then X is a sample from the probability measure [mu]X . Note that the work required to sample X requires only d samples from a univariate induced measure. The procedure above is essentially as described by the authors in [3]; this paper gives a concrete computational method to sample from [[mu].sub.[LAMBDA]] for a relatively general class of measures [mu] (i.e., those formed by arbitrary finite tensor products of Jacobi, half-line Freud, and/or Freud univariate measures). Thus, the algorithms in this paper along with the specifications (5.2) allow one to perform optimal discrete least-squares using Monte Carlo sampling for approximation with multivariate polynomials.

5.2. Weighted equilibrium measures. On [R.sup.d], consider the special case d[mu](x) = exp([-||x||.sup.2]), where || * || is the Euclidean norm on [R.sup.d]. The weighted equilibrium measure [mu]* is a probability measure that is the weak limit of the summations

(5.4) [mathematical expression not reproducible]

where [right arrow] denotes weak (distributional) convergence. The form for [mu]* is not currently known, but the authors in [19] conjecture that [mu]* has support on the unit ball with density

(5.5) [mathematical expression not reproducible]

where [mathematical expression not reproducible]. If X on [R.sup.d] is distributed according to gd, then the cumulative distribution function associated to ||X|| is

(5.6) [mathematical expression not reproducible]

[mathematical expression not reproducible]

where the [r.sup.d-1] factor in the integrand is the [R.sup.d] Jacobian factor for integration in spherical coordinates, and K is the associated normalization constant. Note that the cumulative distribution function Gd is a mapped (normalized) incomplete Beta function with parameters a = d/2 and b = 1 + d/2,

[mathematical expression not reproducible]

where B(*, *) is the Beta function. With d = 1, the veracity of this limit is known [22, Theorem 5.1]. Using the algorithms in this paper, we can empirically test the conjecture.

Precisely, defining [mathematical expression not reproducible], the conjecture for (5.4) reads

[mathematical expression not reproducible]

Our procedure for testing this conjecture is as follows: for a fixed d and large n, we generate M iid samples [mathematical expression not reproducible] distributed according to [mathematical expression not reproducible], and compute the empirical distribution function associated with the ensemble of scalars [mathematical expression not reproducible] We show in Figure 5.1 that indeed for large n that empirical distributions associated with these ensembles match very closely with the distribution function [G.sub.d](r), giving evidence that supports, but does not prove, the conjecture (5.5).

6. Conclusions. We have developed a robust algorithm for the evaluation of induced polynomial distribution functions associated with a relatively wide class of continuous univariate measures. Our algorithms cover all classical orthogonal polynomial measures, and are equally applicable on bounded or unbounded domains. The algorithm leverages several properties of orthogonal polynomials in order to attain stability and accuracy, even for extremely large values of parameters defining the measure or polynomial degree. All computations have been tested up to degree n = 1000 and were found to be stable. The ability to evaluate induced distributions allows the possibility to exactly sample from additive mixtures of these measures. Such additive mixtures define sampling densities that are known to be optimal for multivariate discrete least-squares polynomial approximation algorithms, and allow us to provide supporting empirical evidence for an asymptotic conjecture involving weighted pluripotential equilibrium measures.

Appendix A. Auxiliary recurrences.

For some algorithmic tasks that we consider, the three-term recurrence (2.1) for the pn does not provide a suitable computational procedure due to floating-point under- and over-flow. This happens in two particular cases:

* If x is far outside the convex hull of supp [mu], then [p.sub.n](x) becomes very large and causes numerical overflow (the quantity grows like [x.sup.n]). We will need to evaluate pn(x)/[p.sub.n-1](x) for large x and potentially large n. (When supp[mu] is infinite, one can interpret "far outside supp [mu]" to be defined using the potential-theoretic Mhaskar-Rakhmanov-Saff numbers for a [mathematical expression not reproducible])

* When x is inside the convex hull of supp [mu], we need to evaluate [mathematical expression not reproducible].

For large enough n, direct computation causes numerical overflow.

We emphasize that (2.1) is quite stable and sufficient for most practical computations requiring evaluations of orthogonal polynomials. The situations we describe above (which occur in this paper) are relatively pathological.

A.1. Ratio evaluations. We consider the first case described above. With n fixed, suppose that either x > max[mathematical expression not reproducible], orx < minp [mathematical expression not reproducible]. Then by the interlacing properties of orthogonal polynomial zeros, pj (x) [not equal to] 0 for all j = 0,... n - 1. In this case, the ratio

(A.1) [mathematical expression not reproducible]

is well-defined, with [r.sub.0] := [p.sub.0]. A straightforward manipulation of (2.1) yields

(A.2) [mathematical expression not reproducible]

(A.2) [mathematical expression not reproducible]

The recurrence (A.2) is a more stable way to compute [r.sub.n](x) when x is very large. In practice we can computationally verify that x lies outside the zero set of [p.sub.n-1] with 0(n) effort (e.g., [21, equation (11)] for a crude but general estimate). In the context of this paper, we require linear modifications only on line 8 of Algorithm 1. Here, [mathematical expression not reproducible] whose support is [-1,1], and in the algorithm x e [- 1,1], so the root (3 - x)/(1 + x) has magnitude 1 or greater. Thus, the root is always outside the zero set of [p.sub.n-1] and [r.sub.n] is always well-defined.

A.2. Normalized polynomials. In the second case, consider a different normalization of [p.sub.n]:

(A.3) [mathematical expression not reproducible]

with [mathematical expression not reproducible]. A manipulation of the three-term recurrence relation (2.1) yields the following recurrence for [C.sub.n]:

(A.4a) [mathematical expression not reproducible]

(A.4b) [mathematical expression not reproducible]

(A.4c) [mathematical expression not reproducible]

(A.4d) [mathematical expression not reproducible]

Note that [C.sub.n](x) essentially behaves like rn outside a compact interval containing the zero set of [p.sub.n]; however, [C.sub.n] is well-defined and well-behaved inside this compact interval, unlike [r.sub.n]. The polynomials [p.sub.n] may be reproduced from knowledge of [C.sub.n]:

(A.5) [mathematical expression not reproducible]

To see the veracity of the above, note that

[mathematical expression not reproducible]

and use this to show that [mathematical expression not reproducible] is a telescoping product, verifying that (A.5) is equivalent to (A.3).

Appendix B. Polynomial measure modifications.

We will need to compute recurrence coefficients for the modified measures with densities

[mathematical expression not reproducible]

where we assume that the recurrence coefficients of [mu] are available to us. Here, both [y.sub.0] and [z.sub.0] are some fixed real-valued numbers. In the first case (a linear modification) we assume [y.sub.0] [mathematical expression not reproducible] supp [mu] and choose the sign to ensure that [d.sub.[mu]](x) is positive for x [member of] supp [mu]. Assuming the recurrence coefficients [a.sub.n] and [b.sub.n] for [mu] are known, the problems of computing the recurrence coefficients [a.sub.n] and [b.sub.n] for [mu], and of computing the recurrence coefficients [a.sub.n] and [b.sub.n] for [mu], are well-studied and have constructive computational solutions [9].

We use the auxiliary variables defined in Appendix A to accomplish measure modifications. The linear and quadratic modification recurrence coefficients have the following forms (cf. [18, Section 4]):

(B.1a) [mathematical expression not reproducible]

(B.1b) [mathematical expression not reproducible]

The correction factors for n > 0 are given by

[mathematical expression not reproducible]

For n = 0 they take the special forms

[mathematical expression not reproducible]

Above, [r.sub.n](x) = [r.sub.n] (x;[mu]) and [C.sub.n](x) = [C.sub.n] (x;[mu]) are the functions associated with the measure [mu] and so may be readily evaluated using (A.2) and (A.4).

Note that if we only have a finite number of recurrence coefficients, [mathematical expression not reproducible] for [mu], then a linear modification can only compute modified coefficients up to index M - 1, and a quadratic modification can only compute coefficients up to index M - 2.

Appendix C. Freud and half-line Freud recurrence coefficients. For both cases of Freud measures with [alpha] = 2 (generalized Hermite polynomials), and generalized Freud measures with [alpha] = 1 (generalized Laguerre polynomials), explicit forms for the recurrence coefficients are known. However, the situation is more complicated for other values of [alpha].

We give an extension of Lemma 3.3: Recurrence coefficients of generalized Freud weights may be computed from those of Freud weights.

LEMMA C. 1. Let parameters p[alpha], ([alpha], [rho]) define a Freud weight having recurrence coefficients [mathematical expression not reproducible]. (The [a.sub.n] coefficients vanish because the weight is even.) Define p[[alpha]*,[rho ] *]) and p [alpha]**,[rho ] **) as in (3.19), along with the associated generalized Freud measures [mathematical expression not reproducible]) and [mathematical expression not reproducible]. and their recurrence coefficients [mathematical expression not reproducible] and [mathematical expression not reproducible] respectively. Then, for all n:

(C.1a) [mathematical expression not reproducible]

(C.1b) [mathematical expression not reproducible]


(C.2a) [mathematical expression not reproducible]

(C.2b) [mathematical expression not reproducible]

This result implies that one may either use Freud weight recurrence coefficients to compute half-line Freud weight recurrence coefficients, or vice versa. The proof is a result of Lemma 3.3 along with manipulations of the three-term recurrence relation (2.1).


[1] T. BLOOM, L. BOS, N. LEVENBERG, AND S. WALDRON, On the convergence of optimal measures, Constr.

Approx., 32 (2010), pp. 159-179.

[2] A. COHEN, M. A. DAVENPORT, AND D. LEVIATAN, On the stability and accuracy of least squares approximations, Found. Comput. Math., 13 (2013), pp. 819-834.

[3] A. COHEN AND G. MIGLIORATI, Optimal weighted least-squares methods, SMAI J. Comput. Math., 3 (2017), pp. 181-203.

[4] G. FREUD, Orthogonal Polynomials, Pergamon Press, Oxford, 1971.

[5] W. GAUTSCHI, The interplay between classical analysis and (numerical) linear algebra - A tribute to Gene H. Golub, Electron. Trans. Numer. Anal., 13 (2002), pp. 119-147.

[6] , Orthogonal Polynomials: Computation and Approximation, Oxford University Press, New York, 2004.

[7] , Variable-precision recurrence coefficients for nonstandard orthogonal polynomials, Numer. Algorithms, 52 (2009), pp. 409-418.

[8] W. GAUTSCHI AND S. LI, A set of orthogonal polynomials induced by a given orthogonal polynomial, Aequationes Math., 46 (1993), pp. 174-198.

[9] G. H. GOLUB, Some modified matrix eigenvalue problems, SIAM Rev., 15 (1973), pp. 318-334.

[10] G. H. GOLUB AND J. H. WELSCH, Calculation of Gauss quadrature rules, Math. Comp. 23 (1969), 221-230.

[11] J. HAMPTON AND A. DOOSTAN, Coherence motivated sampling and convergence analysis of least squares polynomial Chaos regression, Comput. Methods Appl. Mech. Engrg., 290 (2015), pp. 73-97.

[12] J. D. JAKEMAN, A. NARAYAN, AND T. ZHOU, A generalized sampling and preconditioning scheme for sparse approximation of polynomial chaos expansions, SIAM J. Sci. Comput., 39 (2017), pp. A1114-A1144.

[13] S. LI, Construction and computation of a new set of orthogonal polynomials, in Applications and Computation of Orthogonal Polynomials (Oberwolfach, 1998), W. Gautschi, G. Opfer, and G. H. Golub, eds., vol. 131 of Internat. Ser. Numer. Math., Birkhauser, Basel, 1999, pp. 145-151.

[14] D. S. LUBINSKY, H. N. MHASKAR, AND E. B. SAFF, A proof of Freud's conjecture for exponential weights, Constr. Approx., 4 (1988), pp. 65-83.

[15] A. P. MAGNUS, Freud's equations for orthogonal polynomials as discrete Painleve equations, in Symmetries and Integrability of Difference Equations (Canterbury, 1996), P. A. Clarkson and F. W. Nijhoff, eds., vol. 255 of London Math. Soc. Lecture Note Ser., Cambridge Univ. Press, Cambridge, 1999, pp. 228-243.

[16] H. N. MHASKAR AND E. B. SAFF, Where does the sup norm of a weighted polynomial live? (A generalization of incomplete polynomials), Constr. Approx., 1 (1985), pp. 71-91. [17] A. NARAYAN, akilnarayan/induced-distributions: Initial beta release, Software, 2017.

[18] A. NARAYAN AND J. S. HESTHAVEN, Computation of connection coefficients and measure modifications for orthogonal polynomials, BIT, 52 (2012), pp. 457-483.

[19] A. NARAYAN, J. D. JAKEMAN, AND T. ZHOU, A Christoffel function weighted least squares algorithm for collocation approximations, Math. Comp., 86 (2017), pp. 1913-1947.

[20] P. G. NEVAI, Orthogonal Polynomials, American Mathematical Society, Providence, 1979.

[21] P. G. NEVAI AND J. S. DEHESA, On asymptotic average properties of zeros of orthogonal polynomials, SIAM J. Math. Anal., 10 (1979), pp. 1184-1192.

[22] E. B. SAFF AND V. TOTIK, Logarithmic Potentials with External Fields, Springer, Berlin, 1997.

[23] B. SIMON, Ratio asymptotics and weak asymptotic measures for orthogonal polynomials on the real line, J. Approx. Theory, 126 (2004), pp. 198-217.

[24] G. SZEGO , Orthogonal Polynomials, 4th ed., American Mathematical Society, Providence, 1975.

[25] W. VAN ASSCHE, Discrete Painleve equations for recurrence coefficients of orthogonal polynomials, in Difference Equations, Special Functions and Orthogonal Polynomials, S. Elaydi, J. Cushing, R. Lasser, V. Papageorgiou, A. Ruffing, and W. Van Assche, eds., World Scientific, New Jersey, 2005, pp. 687-725.

AKIL NARAYAN ([dagger])

Dedicated to Walter Gautschi on the occasion of his 90th birthday

(*) Received April 27, 2017. Accepted September 25, 2018. Published online on December 4, 2018. Recommended by Tom Bella. A. Narayan is partially supported by NSF DMS-1720416, AFOSR FA9550-15-1-0467, and DARPA EQUiPSN660011524053.

([dagger]) Department of Mathematics, and Scientific Computing and Imaging (SCI) Institute, University of Utah, USA (

(1) See Section 2.1 for technical conditions implying this.

(2) Not shown: we have verified this for numerous values of ([alpha] [beta] n).

(3) More formally, it is a potential with "mass" n/[alpha].

DOI: 10.1553/etna_vol50s71
Table 1.1
Classes of measures considered in this paper. [GAMMA](.) is the Euler
Gamma function, and B(.) is the Beta function.

Jacobi           [mathematical expression not reproducible]
Half-line Freud  [mathematical expression not reproducible]
Freud            [mathematical expression not reproducible]

                                       a > - 1
Jacobi           x [member of] [-1,1]  [beta] >-1
                                       a > 0
Half-line Freud  x>0                   [rho]>-1
                                       a > 0
Freud            x[member of]U         [rho]>-1

Jacobi           = [2.sup.[alpha]+[beta]+1]B([beta] + 1,a + 1)
                 [mathematical expression not reproducible]
Half-line Freud
Freud            [mathematical expression not reproducible]
COPYRIGHT 2018 Institute of Computational Mathematics
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2018 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Author:Narayan, Akil
Publication:Electronic Transactions on Numerical Analysis
Article Type:Report
Date:Mar 1, 2018

Terms of use | Privacy policy | Copyright © 2019 Farlex, Inc. | Feedback | For webmasters