Computation of Induced Orthogonal Polynomial Distributions

TL;DR

This paper introduces a stable, spectrally-accurate algorithm for computing distribution functions of induced orthogonal polynomial measures, enabling efficient sampling crucial for advanced polynomial approximation methods.

Contribution

The authors present a novel, robust algorithm for computing and sampling from induced orthogonal polynomial distributions, with implementation available as open-source code.

Findings

Stable computation for polynomial degrees up to 1000

Efficient sampling method for induced distributions

Application to multivariate polynomial approximation

Abstract

We provide a robust and general algorithm for computing distribution functions associated to induced orthogonal polynomial measures. We leverage several tools for orthogonal polynomials to provide a spectrally-accurate method for a broad class of measures, which is stable for polynomial degrees up to at least degree 1000. Paired with other standard tools such as a numerical root-finding algorithm and inverse transform sampling, this provides a methodology for generating random samples from an induced orthogonal polynomial measure. Generating samples from this measure is one ingredient in optimal numerical methods for certain types of multivariate polynomial approximation. For example, sampling from induced distributions for weighted discrete least-squares approximation has recently been shown to yield convergence guarantees with a minimal number of samples. We also provide publicly-available code that implements the algorithms in this paper for sampling from induced distributions.

Figures

Equations

$$\mathrm{d}\tilde{\mu }(x)$$

$$\mathrm{d}\tilde{\tilde{\mu }}(x)$$

$$F_n\left(x\right)$$

$$F_n\left(x_0\right)\approx \frac{1}{2}.$$

$${\mu }_{J,n}^{(\alpha ,\beta )}\left([-1,x]\right)={\mu }_n\left([-1,x]\right)=F_n(x)$$

$$p_n^2(t)\mathrm{d}{\mu }_J^{(\alpha ,\beta )}(t)={\left[p_n(t){\left(\mathrm{d}{\mu }_J^{(\alpha /2n,\beta /2n)}(t)\right)}^n\right]}^2.$$

$$\theta$$

$$x_0\left({\mu }^{(\alpha ,\beta )},n\right)$$

$$A$$

$$F_n(x)$$

$$={\left(\frac{x+1}{2}\right)}^{\beta +1}\frac{c_J^{(0,\beta )}}{c_J^{(\alpha ,\beta )}}{\int }_{-1}^1{\left(2-\frac{1}{2}(u+1)(x+1)\right)}^{\alpha -A}{U}_{2n+A}(u)\mathrm{d}{\mu }^{(0,\beta )}(u),$$

$$U_{2n+A}(u)$$

$$={\gamma }_n^2{\left(\frac{x+1}{2}\right)}^{2n+A}\underset{(a)}{\underbrace{\prod _{k=1}^A\left[\left(\frac{3-x}{1+x}\right)-u\right]}}\underset{(b)}{\underbrace{\prod _{j=1}^n{\left(u-\left(\frac{2}{x+1}\left(x_{j,n}+1\right)-1\right)\right)}^2}},$$

$$\mathrm{d}{\tilde{\mu }}_n(u)=U_{2n+A}(u)\mathrm{d}{\mu }^{(0,\beta )}(u),$$

$$F_n(x)={\left(\frac{x+1}{2}\right)}^{\beta +1}\frac{c_J^{(0,\beta )}}{c_J^{(\alpha ,\beta )}}{\int }_{-1}^1{\left(2-\frac{1}{2}(u+1)(x+1)\right)}^{\alpha -A}\mathrm{d}{\tilde{\mu }}_n(u).$$

$$I(x)$$

$${\hat{F}}_n(x)$$

$${\hat{F}}_n^{(\alpha ,\beta ),c}\left(x\right)={\int }_x^1{\left[p_n^{(\alpha ,\beta )}(t)\right]}^2\mathrm{d}{\mu }^{(\alpha ,\beta )}(t)={\int }_{-1}^{-x}{\left[p_n^{(\beta ,\alpha )}(t)\right]}^2\mathrm{d}{\mu }^{(\beta ,\alpha )}(t)={\hat{F}}_n^{(\beta ,\alpha )}(-x).$$

$$\left|F_n(x)-{\hat{F}}_n(x)\right|\le C(\alpha ,\beta ,n,M)\prod _{j=0}^M{b}_j\left({\tilde{\mu }}_n\right),$$

$$C(\alpha ,\beta ,n,M)=\frac{2^{\beta +1-A}}{(\beta +1)B(\beta +1,\alpha +1)}{\left(\frac{x_0(n)+1}{4}\right)}^{2M+\beta +1}$$

$$\left|F_n(x)-{\hat{F}}_n(x)\right|\le C(\alpha ,\beta ,n,M)\prod _{j=0}^M{b}_j\left({\tilde{\mu }}_n\right)\lesssim C^{\prime \prime }(\alpha ,\beta ,x)4^{-2M},$$

$$\left|{\int }_{-1}^{1}f(u)\mathrm{d}{\tilde{\mu }}_n(u)-\sum _{j=1}^M{\tilde{w}}_{m}f\left({\tilde{u}}_m\right)\right|=\frac{f^{(2M)}(\tau )}{(2M)!}{\int }_{-1}^1\prod _{j=1}^M{\left(u-{\tilde{u}}_m\right)}^2\mathrm{d}{\tilde{\mu }}_n(u)$$

$${\int }_{-1}^1\prod _{j=1}^M{\left(u-{\tilde{u}}_m\right)}^2\mathrm{d}{\tilde{\mu }}_n(u)=\prod _{j=0}^M{b}_j\left({\tilde{\mu }}_n\right){\int }_{-1}^1{p}_M^2\left(u;{\tilde{\mu }}_n\right)\mathrm{d}{\tilde{\mu }}_n(u)=\prod _{j=0}^M{b}_j\left({\tilde{\mu }}_n\right).$$

$$\left|f^{(2M)}(\tau )\right|$$

$$\le {\left(\frac{x_0+1}{2}\right)}^{2M}{2}^{\alpha -A-2M}\prod _{j=0}^{2M-1}|\alpha -A-j|=2^{\alpha -A}{\left(\frac{x_0+1}{4}\right)}^{2M}\prod _{j=0}^{2M-1}|\alpha -A-j|$$

$$\left|F_n(x)-{\hat{F}}_n(x)\right|$$

$$\le {\left(\frac{x+1}{2}\right)}^{\beta +1}\frac{c_J^{(0,\beta )}}{c_J^{(\alpha ,\beta )}}{2}^{\alpha -A}{\left(\frac{x_0+1}{4}\right)}^{2M}\frac{{\prod }_{j=0}^{2M-1}|\alpha -A-j|}{(2M)!}\prod _{j=0}^M{b}_j\left({\tilde{\mu }}_n\right)$$

$$\left|F_n(x)-{\hat{F}}_n(x)\right|$$

$$\mathrm{d}\mu (x)$$

$$\theta$$

$$k_n\left(\alpha \right)=k_n:=n^{1/\alpha }{\left(\frac{2\sqrt{\pi }\Gamma (\alpha )}{\Gamma \left(\alpha +\frac{1}{2}\right)}\right)}^{1/\alpha }.$$

$$p_n(x){\left(\mathrm{d}{\mu }_{HF}^{\left(\alpha ,\rho \right)}(x)\right)}^{1/2}$$

$$\stackrel{u=x^{\alpha }}{=}p\left(u^{1/\alpha }\right){\left[u^{\rho /2n\alpha }\exp \left(-\frac{1}{2n}u\right)\right]}^n$$

$$=q_{n/\alpha }(u){\left[u^{\rho /2n}\exp \left(-\frac{\alpha }{2n}u\right)\right]}^{n/\alpha },$$

Full text

x∈R

Note that $C_{n}(x)/r_{n}(x)\sim 1$ for large enough $|x|$. (Appendix B) Polynomial measure modifications: given $\mu$ and its associated three-term recurrence coefficients $a_{n}$ and $b_{n}$, computation of the three-term recurrence coefficients associated with the measures $\widetilde{\mu}$ and $\widetilde{\widetilde{\mu}}$, defined as

[TABLE]

The $\pm$ sign in $\widetilde{\mu}$ is chosen so that $\widetilde{\mu}$ is a positive measure.

1. Evaluation of $F_{n}$

This section develops computational algorithms for the evaluation of the induced distribution $F_{n}$ defined in (LABEL:eq:mun-def). These algorithms depend fairly heavily on the form of a Lebesgue density $\mathrm{d}\mu(x)$ (i.e., a positive weight function) for the measure $\mu$ on $\mathbbm{R}$, but the ideas can be generalized to various measures. We consider the classes of weights enumerated in Table LABEL:tab:measures:

• •

(Jacobi weights) $\mathrm{d}\mu_{J}(x)=(1-x)^{\alpha}(1+x)^{\beta}$ for $x\in[-1,1]$ with parameters $\alpha,\beta>-1$.

• •

(Freud weights) $\mathrm{d}\mu_{F}(x)=|x|^{\rho}\exp\left(-|x|^{\alpha}\right)$ for $x\in\mathbbm{R}$ with parameters $\alpha>0$, $\rho>-1$.

• •

(half-line Freud weights) $\mathrm{d}\mu_{HF}(x)=x^{\rho}\exp\left(-x^{\alpha}\right)$ for $x\in[0,\infty)$ with parameters $\alpha>0$, $\rho>-1$.

The induced distribution $F_{n}$ 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

[TABLE]

where $\widehat{F}_{n}(x)$ and $\widehat{F}^{c}_{n}(x)$ represent computational approximations to $F_{n}(x)$ and $F^{c}_{n}(x)$, respectively, and are the outputs from the algorithms that we will develop. A pictorial description of this is given in Figure 1 (left). The constant $x_{0}$ is ideally $F_{n}^{-1}(0.5)$; since we cannot know this value a priori, we use potential-theoretic arguments to compute a value $x_{0}=x_{0}\left(\mu,n\right)$ approximating the median of $F_{n}$,

[TABLE]

For our choice of $x_{0}$, we can provide no estimates for this approximate equality, but empirical evidence in Figure 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_{0}$ and detailing algorithms for $\widehat{F}_{n}$ and $\widehat{F}_{n}^{c}$.

1.1. Jacobi weights

We consider computing induced distribution functions for Jacobi measures $\mu^{(\alpha,\beta)}_{J}$ as defined in Table LABEL:tab:measures. When circumstances are clear, we will write $\mu_{J}^{(\alpha,\beta)}=\mu$ to avoid notational clutter. We seek the distribution function of the induced measure,

[TABLE]

To compute $F_{n}$, we specify an approximate median $x_{0}=x_{0}\left(\mu,n\right)$ satisfying (4), and construct algorithmic procedures to evaluate $\widehat{F}_{n}(x)\approx F_{n}(x)$ (for $x\leq x_{0}$) and $\widehat{F}_{n}^{c}(x)\approx F_{n}^{c}(x)$ (for $x>x_{0}$). Having specified these, we use (3) to compute our approximation to $F_{n}$.

1.1.1. Approximating $x_{0}(n)$

With $\alpha$, $\beta>-1$ and $n\in\mathbbm{N}$ fixed, consider the measure $\mu_{J}^{\left(\alpha/2n,\beta/2n\right)}$. Note that this measure is still a Jacobi measure since $\frac{\alpha}{2n}>-1$ and $\frac{\beta}{2n}>-1$. We may rewrite the integrand in (5)

[TABLE]

The quantity under the square brackets on the right-hand side is, in the language of potential theory, a weighed 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_{0}$ of the induced measure $\mu_{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 $\mathrm{d}\mu_{J}^{(\alpha/2n,\beta/2n)}$ on $[-1,1]$. These numbers for this weight are computed explicitly in [23, pp 206-207]. When $\alpha,\beta\geq 0$, this support interval is $\left[\theta-\Delta,\theta+\Delta\right]$, with

[TABLE]

Since this is the interval where most of the “mass” of the integral in (5) lies, we set $x_{0}$ to be the centroid $\theta$ of this interval:

[TABLE]

Note that the definitions of $\theta,\Delta$ require $\alpha$ and $\beta$ to be non-negative. Without mathematical justification, we extend the formula (6) to valid negative values of $\alpha,\beta$ as well. Figure 1 compares $x_{0}$ and $F_{n}^{-1}\left(0.5\right)$ for certain choices of $\alpha$ and $\beta$.

1.1.2. Computing $F_{n}(x)$

First assume that $x\leq x_{0}$, with $x_{0}$ defined in (6), and define $A\in\mathbbm{N}_{0}$ as

[TABLE]

where $\lfloor\cdot\rfloor$ is the floor function. We transform the integral (5) over $[-1,x]$ onto the standard interval $[-1,1]$ via the substitution $u=\frac{2}{x+1}(t+1)-1$:

[TABLE]

where $U_{2n+A}(u)$ is a degree-$(2n+A)$ polynomial given by

[TABLE]

where $\left\{x_{j,n}\right\}_{j=1}^{n}$ are the $n$ zeros of $p_{n}(\cdot)$. Since we explicitly know the polynomial roots of $U_{2n+A}$, we can absorb the term marked $(a)$ into the measure $\mathrm{d}\mu^{(0,\beta)}(u)$ via $A$ linear modifications, and we can likewise absorb the term $(b)$ via $n$ quadratic modifications. (See Appendix B.) Thus, define

[TABLE]

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,

[TABLE]

The integrand above has a root ($\alpha>0$) or singularity ($\alpha<0$) at $u=\frac{3-x}{x+1}=1+2\frac{1-x}{1+x}\geq 1+2\frac{1-x_{0}}{1+x_{0}}$; 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$ $\widetilde{\mu}_{n}$-Gaussian quadrature to efficiently evaluate it. With $\left(\widetilde{u}_{m},\widetilde{w}_{m}\right)_{m=1}^{M}$ denoting the nodes and weights, respectively, of this quadrature rule, we compute

[TABLE]

The entire procedure is summarized in Algorithm 1. In order to compute $\widehat{F}^{c}_{n}(x)$ for $x>x_{0}$, we use symmetry. Since $x\leftarrow-x$ interchanges the parameters $\alpha$ and $\beta$, then

[TABLE]

Note that if $x>x_{0}\left(\mu^{(\alpha,\beta)},n\right)$, then $-x<x_{0}\left(\mu^{(\beta,\alpha)},n\right)$. Thus, $\widehat{F}_{n}^{c}$ can be computed via the same algorithm for $\widehat{F}_{n}$, but with different values for $\alpha$, $\beta$, and $x$. This is also shown in Algorithm 1.

Theorem 1.1**.**

With $\mu^{(\alpha,\beta)}$, $n\in\mathbbm{N}$, and $x\in[-1,1]$ all given, assume that $x\leq x_{0}$ with $x_{0}$ as in (6). Then the output $\widehat{F}_{n}(x)$ from Algorithm 1 using an $M$-point quadrature rule satisfies

[TABLE]

where $b_{j}\left(\widetilde{\mu}_{n}\right)$ are the $b_{j}$ three-term recurrence coefficients associated to the $x$-dependent measure $\widetilde{\mu}_{n}$ defined in (7). The constant $C$ is

[TABLE]

Note that $C(\alpha,\beta,n,M)$ on the right-hand side of (11) is explicitly computable once $\alpha$, $\beta$, $M$, and $n$ are fixed, and only the product involving $b_{j}\left(\widetilde{\mu}_{n}\right)$ depends on $x$. Furthermore, this last quantity is explicitly computed in Algorithm 1, so that a rigorous error estimate for the algorithm can be computed before its termination. Since $\frac{x_{0}+1}{4}<\frac{1}{2}$, then the estimate (11) also hints at exponential convergence of the quadrature rule strategy, assuming that the factors $b_{j}\left(\widetilde{\mu}_{n}\right)$ can be bounded or controlled. We cannot provide these bounds, although we do know the asymptotic behavior $b_{j}\left(\widetilde{\mu}_{n}\right)\rightarrow\frac{1}{4}$ as $j\rightarrow\infty$ [21, Remark 3.1.10]. Also, since $x_{0}\in[-1,1]$ for any $n>0$ then, uniformly in $n>0$, $C(\alpha,\beta,n,M)\leq C^{\prime}(\alpha,\beta)4^{-M}$. Thus, for a fixed $x$ we expect that the estimate in Theorem 1.1 behaves like

[TABLE]

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

Proof.

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: $\widetilde{u}_{m}$ and $\widetilde{w}_{m}$ denote the $M$-point $\widetilde{\mu}_{n}$-Gaussian quadrature nodes and weights, respectively. We start with the Corollary to Theorem 1.48 in [7], stating that if $f(\cdot)$ is infinitely differentiable on $[-1,1]$, then

[TABLE]

for some $\tau\in(-1,1)$. Noting that since $\widetilde{u}_{m}$ are the zeros of the degree-$M$ $\widetilde{\mu}_{n}$-orthogonal polynomial, then

[TABLE]

From (8), the integral we wish to approximate has integrand $f(u)=\left(2-\frac{1}{2}(u+1)(x+1)\right)^{\alpha-A}$. Then for any $\tau\in(-1,1)$ and any $x\leq x_{0}$,

[TABLE]

Then

[TABLE]

Since $\alpha-A\in(-1,1)$, then

[TABLE]

The result follows by direct computation of $2^{\alpha+\beta+1-A}\frac{c_{J}^{(0,\beta)}}{c_{J}^{(\alpha,\beta)}}$. ∎

We verify spectral convergence of the scheme empirically in Figure 2, 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 $\alpha$ or $\beta$. Based on these results, taking $M=10$ appears sufficient uniformly over all $(\alpha,\beta,n)$.111Not shown: we have verified this for numerous values of $(\alpha,\beta,n)$. Finally, extending the estimate in (11) to the case $x>x_{0}$ can be accomplished by permuting $\alpha$ and $\beta$ as is done for that case in Algorithm 1.

1.2. Half-line Freud weights

In this section we consider the half-line Freud measure $\mu^{(\alpha,\rho)}_{HF}$ as defined in Table LABEL:tab:measures. These algorithms require the recurrence coefficients for $\mu_{HF}$; these coefficients are in general not easy to compute when $\alpha\neq 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 [16, 1]. In our computations we use the methodology of [8] to compute Freud weight recurrence coefficients (and hence half-line Freud coefficients). We note that the methodology of [8] is computationally onerous: For a fixed $\alpha$ and $\rho$, it required a day-long computation to obtain 500 recurrence coefficients. Again in this subsection we write $\mu^{(\alpha,\rho)}_{HF}=\mu$ and similarly for $\mu_{n}$, $F$, $F_{n}$, etc. We restore these super- and subscripts when ambiguity arises without them. We accomplish computation of $F_{n}$ 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 $\textrm{supp}\mu$ (which is $x=0$ here), or to the right-hand end of the $\textrm{supp}\mu$ (which is $x=\infty$). We determine this delineation again by means of potential theory.

1.2.1. Computation of $x_{0}$

As with the Jacobi case, we take $x_{0}$ to be the midpoint of the “essential” support for $p_{n}^{2}(x)\mathrm{d}\mu_{HF}^{(\alpha,\rho)}$, the latter of which is approximately the support of the weighted equilibrium measure associated to $\left[\mathrm{d}\mu_{HF}^{\left(\alpha,\rho\right)}\right]^{1/2n}$. 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

[TABLE]

is the interval $[\theta-\Delta,\theta+\Delta]$, with these values given by [23]

[TABLE]

This interval is the “essential” support for any function of the form $p_{n}(x)\left(\mathrm{d}\mu(x)\right)^{n}$ where $p_{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 $\sqrt{\mathrm{d}\mu_{HF}^{(\alpha,0)}}$ in this case is the interval $[0,k_{n}(\alpha)]$, where

[TABLE]

are the Mhaskar-Rakhmanov-Saff numbers for $\sqrt{\mu_{HF}^{(\alpha,0)}}$ [17].

We now derive our approximation for the case of general $\alpha$, $\rho$. To approximate where $p_{n}\sqrt{\mathrm{d}\mu_{HF}^{(\alpha,\rho)}}$ is supported, consider

[TABLE]

where we have introduced $q_{n/\alpha}$, which is a “polynomial” of “degree” $n/\alpha$222More formally, it is a potential with “mass” $n/\alpha$.. Note that in the variable $u$, the weight function under square brackets in the last expression is of the form (12). Concepts in potential theory extend to generalized notions of polynomial degree, and so we may apply our formulas for $\theta$ and $\Delta$ with $s=\frac{\rho}{2n}$ and $\lambda=\frac{\alpha}{2n}$. These formulas imply that the “essential” support for the $u$ variable is

[TABLE]

Therefore, to obtain appropriate limits on the variable $x$, we raise the endpoints $\theta\pm\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$

[TABLE]

and compare this with the exact value $k_{n}(\alpha)$ computed above. We note that while $k_{n}(\alpha)\sim n^{1/\alpha}$ matches the $n$-behavior of (13), the constant is wrong. We thus multiply the endpoints $\left(\theta\pm\Delta\right)^{1/\alpha}$ by the appropriate constant to match the $\rho=0$ behavior of $k_{n}(\alpha)$. The net result then, for arbitrary $\alpha,\rho$, is the approximation

[TABLE]

Figure 3 compares the intervals demarcated by $a_{-}$ and $a_{+}$ versus $F_{n}^{-1}\left([0.01,0.99]\right)$, the latter of which contains “most” of the support of $F_{n}$.

1.2.2. Computing $\widehat{F}_{n}(x)$

First assume that $x\leq x_{0}$. Then

[TABLE]

We recognize a portion of the integrand as a Jacobi measure, and use successive measure modifications to define $\widetilde{\mu}_{n}$:

[TABLE]

The recurrence coefficients of $\widetilde{\mu}_{n}$ can be computed via polynomial measure modifications on the roots $u_{j,n}=\frac{2x_{j,n}}{x}-1$, where $x_{j,n}$ are the roots of $p_{n}\left(\cdot\right)$. A detailed algorithm is given in Algorithm 2.

Now assume that $x\geq x_{0}$. We compute $F_{n}^{c}$ directly in a similar fashion as we did for $F_{n}$. We have

[TABLE]

We again use this to define a new measure $\widetilde{\mu}_{n}$ and an associated $M$-point Gauss quadrature $\left(u_{m},w_{m}\right)_{m=1}^{M}$. The recurrence coefficients for $\widetilde{\mu}_{n}$ are computable via polynomial measure modifications. This results in the approximation

[TABLE]

A more detailed algorithm is given in Algorithm 3. Of course, once $F_{n}^{c}$ is computed we may compute $F_{n}(x)=1-F_{n}^{c}(x)$.

Errors between $F_{n}$ and computational approximations $\widehat{F}_{n}$ are shown in Figure 4. 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 4. 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 $x_{0}(n)\equiv 50$. Whereas for $\alpha\neq 1$ our tests suggest that $M=n+10$ is sufficient to achieve good accuracy, and we take $x_{0}$ as the average of $a_{\pm}$ as given in (14).

1.3. Freud weights

Finally, consider the Freud measure $\mu_{F}^{(\alpha,\rho)}$ defined in Table LABEL:tab:measures. 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 [15]. 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 1.2**.**

Let parameters $(\alpha,\rho)$ define a Freud weight and associated measure $\mu_{F}^{(\alpha,\rho)}$. Then for $x\leq 0$,

[TABLE]

For $x\geq 0$, we have

[TABLE]

Note that with expressions (21) and (22), an algorithm for computing $F_{n}\left(\cdot;\mu_{F}\right)$ is straightforward to devise utilizing Algorithm 3 for $F_{n}^{c}\left(\cdot;\mu_{HF}\right)$. The result (22) follows easily from the fact that the integrand in (LABEL:eq:mun-def) defining $F_{n}$ is an even function. To prove the main portion of the theorem, expression (21), we require the following result relating Freud orthonormal polynomials to half-line Freud orthonormal polynomials.

Lemma 1.1**.**

Let $\rho>-1$ and $\alpha>1$ be parameters that define a Freud measure $\mu^{(\alpha,\rho)}_{F}$ with associated orthonormal polynomial family $p_{n}(x)=p_{n}\left(x;\mu_{F}^{(\alpha,\rho)}\right)$. Define two sets of half-line Freud parameters $\left(\alpha_{\ast},\rho_{\ast}\right)$ and $\left(\alpha_{\ast\ast},\rho_{\ast\ast}\right)$ and the corresponding half-line Freud measures and polynomials:

[TABLE]

Also define the constant

[TABLE]

Then, for all $n\geq 0$,

[TABLE]

Proof.

The following equalities may be verified via direct computation using the definitions (23) and (24) along with the expressions in Table LABEL:tab:measures:

[TABLE]

The proof of this lemma relies on the change of measure $t\mapsto x^{2}$. We have

[TABLE]

This relation shows that the family $\left\{p_{\ast,n}\left(x^{2}\right)\right\}_{n=0}^{\infty}$ are polynomials of degree $2n$ that are orthonormal under a Freud weight with parameters $\alpha=2\alpha_{\ast}$ and $\rho=2\rho_{\ast}+1$. Using nearly the same arguments, but with the family $p_{\ast\ast,n}$, yields the relation

[TABLE]

This relation shows that the family $\left\{hxp_{\ast\ast,n}\left(x^{2}\right)\right\}_{n=0}^{\infty}$ are polynomials of degree $2n+1$ that are orthonormal under the same Freud having parameters $\alpha=2\alpha_{\ast\ast}$ and $\rho=2\rho_{\ast\ast}-1$. We also have that $xp_{\ast\ast,n}\left(x^{2}\right)$ is orthogonal to $p_{\ast,m}\left(x^{2}\right)$ under a(ny) Freud weight for any $n,m$ because of even-odd symmetry. Thus, define

[TABLE]

Then $\left\{P_{n}\right\}_{n=0}^{\infty}$ is a family of degree-$n$ polynomials (with positive leading coefficient) orthonormal under a $(\alpha,\rho)$ Freud weight. Therefore, $P_{n}\equiv p_{n}$. ∎

We can now give the

Proof of Thereom 1.2.

Assume $x\leq 0$. Then

[TABLE]

where $(\alpha_{\ast},\rho_{\ast})$ are as defined in (23a). By Lemma 1.1,

[TABLE]

where $(\alpha_{\ast\ast},\rho_{\ast\ast})$ are defined in (23b). Then if $n$ is even, we have

[TABLE]

where we recall the equalities (25) related $c_{F}$ to $c_{F\ast}$. Similarly, if $n$ is odd we have

[TABLE]

The combination of these results proves (21). ∎

2. Inverting induced distributions

We have discussed at length in previous sections algorithms for computing $F_{n}(x)$ defined in (LABEL:eq:mun-def) for various Lebesgue-continuous measures $\mu$. The central application of these algorithms we investigate in this paper is actually in the evaluation of $F_{n}^{-1}(u)$ for $u\in[0,1]$. We accomplish this by solving for $x$ in the equation

[TABLE]

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

2.1. Computing an initial interval

We use $s_{\pm}$ to denote the (possibly infinite) endpoints of the support of $\mu$:

[TABLE]

Now let $u\in[0,1]$. Our first step in finding $F_{n}^{-1}(u)$ is to compute two values $x_{-}$ and $x_{+}$ such that

[TABLE]

Our procedure for identifying an initial interval containing $F_{n}^{-1}(u)$ leverages the Markov-Stiltjies 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 2.1** (Markov-Stiltjies Inequalties, [25]).**

Let $\mu$ be a probability measure on $\mathbbm{R}$ with an associated orthogonal polynomial family. For any $N\in\mathbbm{N}$, let $\left\{x_{k,N},w_{k,N}\right\}_{k=1}^{N}$ denote the $N$ $\mu$-Gaussian quadrature nodes and weights, respectively. Then:

[TABLE]

Let $\left\{p_{j,n}\right\}_{j=0}^{\infty}$ denote the sequence of polynomial orthonormal under the induced measure $\mu_{n}$. Given $N\in\mathbbm{N}$, let $\left\{z_{k,N},v_{k,N}\right\}_{k=1}^{N}$ denote the $N$-point $\mu_{n}$-Gaussian quadrature rule, i.e., $z_{k,N}$ are the $N$ ordered zeros of $p_{N,n}(x)$, with $v_{k,N}$ the associated weights. Since $\mu_{n}$ is a probability measure, then $\sum_{k=1}^{N}v_{k,N}=1$. As such, given $u\in[0,1]$ we can always find some $m\in\left\{1,\ldots,N\right\}$ such that

[TABLE]

Then, defining $z_{0,N}\equiv s_{-}$ and $z_{N+1,N}\equiv s_{+}$ for all $N$ and $n$, we have

[TABLE]

Since $F_{n}$ is non-decreasing, this is equivalently,

[TABLE]

Thus, if we find an $m$ such that (28) holds, then (27) holds with

[TABLE]

When $\mathrm{supp}\;\mu$ is bounded, the $N$-asymptotic density of orthogonal polynomial zeros on $\mathrm{supp}\;\mu$ guarantees that we can find a bounding interval with endpoints $x_{\pm}$ of arbitrarily small width by taking $N$ sufficiently large. The difficulty is that we therefore require the zeros $z_{k,N}$ and the quadrature weights $v_{k,N}$ of the induced measure, which in turn require knowledge of the three-term recurrence coefficients associated to $\mu_{n}$. These can be easily computed from the coefficients associated to $\mu$; since

[TABLE]

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

2.2. Bisection

For simplicity, the root-finding method we employ to solve (26) is the bisection approach. More sophisticated methods may be used, with the caveat that the derivative of the function, $F^{\prime}(x)=p_{n}^{2}(x)\mathrm{d}\mu(x)$, vanishes wherever $p_{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 (26) starts with an initial guess for an interval $[x_{-},x_{+}]$ containing the root $x$, and iteratively updates this interval via

[TABLE]

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_{\pm}$ lessens the number of evaluations of $F_{n}$ in a bisection approach and thus accelerates overall evaluation of $F_{n}^{-1}$. The overall algorithm for solving (26) is to (i) compute the recurrence coefficients associated with $\mu_{n}$ in (30) via quadratic measure modifications, (ii) compute order-$N$ $\mu_{n}$-Gaussian quadrature nodes and weights $z_{j,N}$ and $v_{j,N}$, respectively, (iii) identify $m$ such that (28) holds so that $x_{\pm}$ may be computed in (29), and (iv) iteratively apply the bisection algorithm with the initial interval defined by $x_{\pm}$ using the evaluation procedures for $F_{n}$ outlined in Section 1.

3. 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\in\mathbbm{N}_{0}^{d}$ denotes a multi-index with components $\lambda=\left(\lambda_{1},\ldots,\lambda_{d}\right)$ and magnitude $|\lambda|=\sum_{j=1}^{d}\lambda_{j}$. A point $x\in\mathbbm{R}^{d}$ has components $x=\left(x_{1},\ldots,x_{d}\right)$, and $x^{\lambda}=\prod_{j=1}^{d}x_{j}^{\lambda_{j}}$. A collection of multi-indices will be denoted $\Lambda\subset\mathbbm{N}_{0}^{d}$; we will assume that $N=\left|\Lambda\right|$ is finite. Let $\mu$ be a tensorial measure on $\mathbbm{R}^{d}$ such that each of its $d$ marginal univariate measures $\mu^{(j)}$, $j=1,\ldots,d$ admits a $\mu^{(j)}$-orthonormal polynomial family $\left\{p_{j,n}\right\}_{n=0}^{\infty}$ on $\mathbbm{R}$, satisfying

[TABLE]

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

[TABLE]

These polynomials are an $L^{2}_{\mathrm{d}\mu}$-orthonormal basis for the subspace $P_{\Lambda}$, defined as

[TABLE]

Under the additional assumption that the index set $\Lambda$ is downward-closed, then $P_{\Lambda}=\mathrm{span}\left\{x^{\lambda}\;\;|\;\;\lambda\in\Lambda\right\}$. We extend our definition of induced polynomials to this tensorial multivariate situation. For any $\lambda\in\Lambda$, the order-$\lambda$ induced measure $\mu_{\lambda}$ is defined as

[TABLE]

where $\mathrm{d}\mu^{(j)}_{\lambda_{j}}$ is the (univariate) order-$\lambda_{j}$ induced measure for $\mu^{(j)}$ according to the definition (LABEL:eq:mun-def). Thus, $\mu_{\lambda}$ is also a tensorial measure.

3.1. Optimal polynomial discrete least-squares

The goal of this section is description of 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 [4] and is based on the foundational matrix concentration estimates for least-squares derived in [3]. Let $f:\mathbbm{R}^{d}\rightarrow\mathbbm{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 $P_{\Lambda}$, we are interested in approximating the $L^{2}_{\mathrm{d}\mu}$-orthogonal projection of $f$ onto $P_{\Lambda}$. This projection is given explicitly by

[TABLE]

One way to approximate the integral defining the coefficients $c^{\ast}_{\lambda}$ is via a Monte Carlo least-squares procedure using $M$ collocation samples of the function $f(x)$. Let $\left\{X_{m}\right\}_{m=1}^{M}$ denote a collection of $M$ independent and identically distributed random variables on $\mathbbm{R}^{d}$, where we leave the distribution of $X_{m}$ unspecified for the moment. A weighted discrete least-squares recovery procedure approximates $c^{\ast}_{\lambda}$ with $c_{\lambda}$, computed as

[TABLE]

where $w_{m}$ are positive weights. One supposes that if the distribution of $X_{m}$ and the weights $w_{m}$ are chosen intelligently, then it is possible to recover the $N$ coefficients $c_{\lambda}$ with a relatively small number of samples $M$; ideally, $M$ should be close to $N$. The analysis in [3] codifies conditions on a required sample count $M$ so that the minimization procedure above is stable, and so that the recovered coefficients $c_{\lambda}$ are “close” to $c_{\lambda}^{\ast}$; these conditions depend on the distribution of $X_{m}$, on $w_{m}$, on $\mu$, and on $P_{\Lambda}$. Since $\mu$ and $P_{\Lambda}$ are specified, the goal here is identification of an appropriate distribution for $X_{m}$ and weight $w_{m}$. Using ideas proposed in [20, 12] the results in [4] show that, in the context of the analysis in [3], the optimal choice of probability measure $\mu_{X}$ for sampling $X_{m}$ and weights $w_{m}$ that achieves a minimal sample count $M$ are

[TABLE]

The precise quantification of the sample count and error estimates can be formulated using an algebraic characterization of (31). Define matrices $\boldsymbol{V}\in\mathbbm{R}^{M\times N}$ and $\boldsymbol{W}\in\mathbbm{R}^{M\times M}$, and vectors $\boldsymbol{c}\in\mathbbm{R}^{N}$ and $\boldsymbol{f}\in\mathbbm{R}^{M}$ as follows:

[TABLE]

where $\lambda(1),\ldots,\lambda(N)$ represents any enumeration of elements in $\Lambda$. We use $\|\cdot\|$ on matrices to denote the induced $\ell^{2}$ norm. The algebraic version of (31) is then to compute $\boldsymbol{c}$ that minimizes the the least-squares residual of $\sqrt{\boldsymbol{W}}\boldsymbol{V}\boldsymbol{c}=\sqrt{\boldsymbol{W}}\boldsymbol{f}$. The following result holds.

Theorem 3.1** ([3, 4]).**

Let $0<\delta<1$, and $r>0$ be given, and define $c_{\delta}\coloneqq\delta+(1-\delta)\log(1-\delta)\in(0,1)$. Draw $M$ iid samples $\left\{X_{m}\right\}_{m=1}^{M}$ from $\mu_{X}$, and let the coefficients $c_{\lambda}$ be those recovered from (31). If

[TABLE]

Then

[TABLE]

The free parameter $r$ is a tunable oversampling rate; $\delta$ represents the guaranteeable proximity of $\boldsymbol{V}^{T}\boldsymbol{W}\boldsymbol{V}$ to $\boldsymbol{I}$. We emphasize that by choosing $\mu_{X}=\mu_{\Lambda}$ with the weights defined as in $\boldsymbol{W}$, then the size of $M$ as required by (33) depends only on the the cardinality $N$ of $\Lambda$, and not on its shape. Furthermore, the criterion $M/\log M\gtrsim 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_{X}$ is an additive mixture of induced measures and can be easily sampled, assuming $\mu_{\lambda}$ can be sampled. Sampling from $\mu_{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 $U_{j}$, $j=1,\ldots,d$ each on the interval $[0,1]$, (3) compute $X\in\mathbbm{R}^{d}$ as

[TABLE]

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 [4]; this paper gives a concrete computational method to sample from $\mu_{\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 (32) allow one to perform optimal discrete least-squares using Monte Carlo sampling for approximation with multivariate polynomials.

3.2. Weighted equilibrium measures

On $\mathbbm{R}^{d}$, consider the special case $\mathrm{d}\mu(x)=\exp(-\|x\|^{2})$, where $\|\cdot\|$ is the Euclidean norm on $\mathbbm{R}^{d}$. The weighted equilibrium measure $\mu^{\ast}$ is a probability measure that is the weak limit of the summations

[TABLE]

The form for $\mu^{\ast}$ is not currently known, but the authors in [20] conjecture that $\mu^{\ast}$ has support on the unit ball with density

[TABLE]

where $C_{d}=\left(\pi\right)^{-(d+1)/2}\Gamma\left(\frac{d+1}{2}\right)$. If $X$ on $\mathbbm{R}^{d}$ is distributed according to $g_{d}$, then the cumulative distribution function associated to $\left\|X\right\|$ is

[TABLE]

where the $r^{d-1}$ factor in the integrand is the $\mathbbm{R}^{d}$ Jacobian factor for integration in spherical coordinates, and $K$ is the associated normalization constant. Note that the cumulative distribution function $G_{d}$ is a mapped (normalized) incomplete Beta function with parameters $a=d/2$ and $b=1+d/2$,

[TABLE]

where $B(\cdot,\cdot)$ is the Beta function. With $d=1$, the veracity of this limit is known . Using the algorithms in this paper, we can empirically test the conjecture. Precisely, defining $\Lambda_{n}\coloneqq\left\{\lambda\in\mathbbm{N}_{0}^{d}\;|\;|\lambda|\leq n\right\}$, then the conjecture for (34) reads

[TABLE]

Our procedure for testing this conjecture is as follows: for a fixed $d$ and large $n$, we generate $M$ iid samples $\left\{X_{m}\right\}_{m=1}^{M}$ distributed according to $\mu_{\Lambda_{n}}$, and compute the empirical distribution function associated with the ensemble of scalars $\left\{\left\|X_{m}\right\|/\sqrt{2n}\right\}_{m=1}^{M}$. We show in Figure 5 that indeed for large $n$ that empirical distributions associated with these ensembles match very closely with the distribution function $G_{d}(r)$, giving evidence that supports, but does not prove, the conjecture (35).

4. 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 (LABEL:eq:three-term-recurrence) for the $p_{n}$ 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 $\textrm{supp}{}\mu$, then $p_{n}(x)$ becomes very large and causes numerical overflow (the quantity grows like $x^{n}$). We will need to evaluate $p_{n}(x)/p_{n-1}(x)$ for large $x$ and potentially large $n$. (When $\textrm{supp}{}\mu$ is infinite, one can interpret “far outside $\textrm{supp}{}\mu$” to be defined using the potential-theoretic Mhaskar-Rakhmanov-Saff numbers for $\sqrt{\mathrm{d}\mu(x)}$.)

• •

When $x$ is inside $\textrm{supp}\mu$, we will need to evaluate $p_{n}^{2}(x)/\sum_{j=0}^{n-1}p_{j}^{2}(x)$. For large enough $n$, a direct computation causes numerical overflow.

We emphasize that (LABEL:eq:three-term-recurrence) 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 p_{n-1}^{-1}(0)$, or $x<\min p_{n-1}^{-1}(0)$. Then by the interlacing properties of orthogonal polynomial zeros, $p_{j}(x)\neq 0$ for all $j=0,\ldots n-1$. In this case, the ratio

[TABLE]

is well-defined, with $r_{0}\coloneqq p_{0}$. A straightforward manipulation of (LABEL:eq:three-term-recurrence) yields

[TABLE]

The recurrence (38) is a more stable way to compute $r_{j}(x)$ when $x$ is very large. In practice we can computationally verify that $x$ lies outside the zero set of $p_{n-1}$ with $\mathcal{O}(n)$ effort (e.g., [22, equation (11)] for a crude but general estimate). In the context of this paper, this condition is always satisfied whenever we require an evaluation of $r_{n}(x)$.

A.2. Normalized polynomials

In the second case, consider a different normalization of $p_{n}$:

[TABLE]

with $C_{0}\equiv p_{0}=\sqrt{1/b_{0}}$. A manipulation of the three-term recurrence relation (LABEL:eq:three-term-recurrence) yields the following recurrence for $C_{n}$:

[TABLE]

Note that $C_{n}(x)$ essentially behaves like $r_{n}$ outside a compact interval containing the zero set of $p_{n}$; however, $C_{n}$ is well-defined and well-behaved inside this compact interval, unlike $r_{n}$. The polynomials $p_{n}$ may be reproduced from knowledge of $C_{n}$:

[TABLE]

Appendix B Polynomial measure modifications

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

[TABLE]

where we assume that the recurrence coefficients of $\mu$ are available to us. Here, both $y_{0}$ and $z_{0}$ are some fixed real-valued numbers. In the first case (a linear modification) we assume $y_{0}\not\in\textrm{supp}\mu$ and choose the sign to ensure that $\mathrm{d}\widetilde{\mu}(x)$ is positive for $x\in\textrm{supp}\mu$. Assuming the recurrence coefficients $a_{n}$ and $b_{n}$ for $\mu$ are known, the problems of computing the recurrence coefficients $\widetilde{a}_{n}$ and $\widetilde{b}_{n}$ for $\widetilde{\mu}$, and of computing the recurrence coefficients $\widetilde{\widetilde{a}}_{n}$ and $\widetilde{\widetilde{b}}_{n}$ for $\widetilde{\widetilde{\mu}}$, are well-studied and have constructive computational solutions [10]. 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. [19, Section 4]):

[TABLE]

The correction factors for $n>0$ are given by

[TABLE]

For $n=0$ they take the special forms

[TABLE]

Above, $r_{n}(x)=r_{n}\left(x;\mu\right)$ and $C_{n}(x)=C_{n}\left(x;\mu\right)$ are the functions associated with the measure $\mu$ and so may be readily evaluated using (38) and (40). Note that if we only have a finite number of recurrence coefficients, $\left\{a_{n},b_{n}\right\}_{n=0}^{M}$ 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 1.1: Recurrence coefficients of generalized Freud weights may be computed from those of Freud weights.

Lemma C.1**.**

Let parameters $(\alpha,\rho)$ define a Freud weight having recurrence coefficients $\left\{b_{n}\right\}_{n=0}^{\infty}$. (The $a_{n}$ coefficients vanish because the weight is even.) Define $(\alpha_{\ast},\rho_{\ast})$ and $(\alpha_{\ast\ast},\rho_{\ast\ast})$ as in (23), along with the associated generalized Freud measures $\mu_{F\ast}^{(\alpha_{\ast},\rho_{\ast})}$ and $\mu_{F\ast}^{(\alpha_{\ast\ast},\rho_{\ast\ast})}$ and their recurrence coefficients $\left\{(a_{\ast,n},b_{\ast,n})\right\}_{n=0}^{\infty}$ and $\left\{(a_{\ast\ast,n},b_{\ast\ast,n})\right\}_{n=0}^{\infty}$, respectively. Then, for all $n$:

[TABLE]

Furthermore,

[TABLE]

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 1.1 along with manipulations of the three-term recurrence relation (LABEL:eq:three-term-recurrence).