Construction and implementation of asymptotic expansions for Laguerre–type orthogonal polynomials

Construction and implementation of asymptotic expansions for Laguerre–type orthogonal polynomials

Daan Huybrechs
Department of Computer Science
KU Leuven, Belgium
   Peter Opsomer (corresponding author)
Department of Computer Science
KU Leuven, Belgium

Laguerre and Laguerre-type polynomials are orthogonal polynomials on the interval with respect to a weight function of the form

The classical Laguerre polynomials correspond to . The computation of higher-order terms of the asymptotic expansions of these polynomials for large degree becomes quite complicated, and a full description seems to be lacking in literature. However, this information is implicitly available in the work of Vanlessen [28], based on a non-linear steepest descent analysis of an associated so-called Riemann–Hilbert problem. We will extend this work and show how to efficiently compute an arbitrary number of higher-order terms in the asymptotic expansions of Laguerre and Laguerre-type polynomials. This effort is similar to the case of Jacobi and Jacobi-type polynomials in a previous paper. We supply an implementation with explicit expansions in four different regions of the complex plane. These expansions can also be extended to Hermite-type weights of the form on , and to general non-polynomial functions using contour integrals. The expansions may be used, e.g., to compute Gauss-Laguerre quadrature rules in a lower computational complexity than based on the recurrence relation, and with improved accuracy for large degree. They are also of interest in random matrix theory.

1 Introduction

We determine asymptotic approximations as of the orthonormal polynomials on with positive leading coefficient, with the weight function

The classical Laguerre polynomials corresponds to , but we aim to provide formulas and results for general functions . The choice corresponds to so-called Freud-type polynomials [17].

The procedure in [28] gives four types of asymptotic expansions: (\Romanbar1) inner asymptotics for near the bulk of the zeros but away from the extreme zeros, (\Romanbar2) outer asymptotics valid for away from the zeros of , (\Romanbar3) boundary asymptotics near the so-called soft edge valid for near the largest zeros and (\Romanbar4) boundary asymptotics near the so-called hard edge for near . We also provide asymptotic expansions for associated quantities such as leading term coefficients as well as recurrence coefficients and of the three term recurrence relation


The methodology of [28] is based on the non-linear steepest descent method by Deift and Zhou [9] for a Riemann–Hilbert problem that is generically associated with orthogonal polynomials by Fokas, Its and Kitaev [10]. This is further detailed in § 2.1. The general strategy is to apply a sequence of transformations , such that the final matrix-valued function is asymptotically close to the identity matrix as or tends to . The asymptotic result for , and subsequently for the polynomials, is obtained by inverting these transformations. The transformations involve a normalization of behaviour at infinity, the so-called ‘opening of a lens’ around the interval of orthogonality, and the introduction of local parametrices in disks around special points like endpoints, which are matched to global parametrices elsewhere in the complex plane. These transformations split into different regions, where different formulas for the asymptotics are valid.

In our case of Laguerre-type polynomials, one also first needs an -dependent rescaling of the axis using the so-called MRS numbers (defined further on in this paper). After this step the roots of the rescaled polynomials accumulate in a fixed and finite interval. We provide an algorithm to obtain an arbitrary number of terms in the expansions, where we set up series expansions using many convolutions that follow the chain of transformations and their inverses in this steepest descent method. While doing this, we keep computational efficiency in mind, as well as the use of the correct branch cuts in the complex plane.

The strategy outlined above and in § 2.1 was also followed in our earlier article about asymptotic expansions of Jacobi–type polynomials [6], which was based on the mathematical analysis of Kuijlaars et al in [16]. Here, we base our results on the analysis in the work of Vanlessen [28]. The main differences in this paper compared to [6] are the following:

  • The analysis of Laguerre-type polynomials on the halfline has an extra step that involves a rescaling via the MRS numbers (see § 3). This leads to fractional powers of in the expansions.

  • There is also a new behaviour near the largest zero (often referred to as a soft edge), captured by the Airy function. This leads to higher order poles in the derivations and thus also to longer formulas for the higher order terms in the expansions. The behaviour near the hard edge at involves Bessel functions, like in the Jacobi case near the endpoints .

  • We obtain more explicit results for polynomials . We will frequently distinguish in this paper between three cases: monomial , general polynomial and a more general analytic function .

Asymptotic expansions can be useful in computations for several reasons. The use of the recurrence relation (1.1) in applications involving Laguerre polynomials results in accumulating roundoff errors and a computation time that is linear in . In contrast, asymptotic expansions become increasingly accurate as the degree becomes large, and the computing time is essentially independent of . Another motivation is the partition function for Laguerre ensembles of random matrices as mentioned in [28] and studied in [29]. This gives the eigenvalue distribution of products of random matrices of a certain type, which could for example arise in stochastic processes (Markov chains) or quantum mechanics.

In this reference, some leading order terms are given explicitly, and we detail the derivation of higher-order terms. The analysis only requires elementary numerical techniques: in particular, there is no need for the evaluation of special functions besides the Airy and Bessel functions. The formulas are implemented in Sage and Matlab and are available on the software web page of our research group [21].

The expressions are also implemented in the Chebfun package in Matlab for computing with functions (see [24, lagpts.m]) and into the Julia package FastGaussQuadrature [25, gausslaguerre.jl], in both cases for the purpose of constructing Gauss-Laguerre quadrature rules with a high number of points in linear complexity. This approach is comparable to a number of modern numerical methods for other types of Gaussian quadrature, that are often based on asymptotics and that lead to a linear complexity as well [12, 2, 14, 1]. A recent paper [4] achieves competitive performance in a general way via nonoscillatory phase functions. Potential improvements to our code in light of this result and a more thorough discussion of the contributions to the computation of Gaussian quadrature rules are future research topics. We do remark here that the construction of Gaussian quadrature rules requires the derivative of the associated orthonormal Laguerre polynomial with positive leading coefficient. This can be obtained from our expansions and the identity , derived from [19, 18.9.23]. Although technical, the expansions can readily be differentiated for general . This paper affirmatively answers the questions raised in the conclusions of [26], namely whether the RH approach can be applied to the fast computation of quadrature rules with generalized Laguerre weights (for the Jacobi case, see [6]), and whether higher order asymptotic expansions can be computed effectively.

As mentioned, the standard Laguerre polynomials correspond to , or equivalently , , and . For this case, asymptotic expansions are given in [5] with explicit expressions for the first terms. We refer the reader to [23, 18] and references therein for more results on asymptotics for the standard Laguerre polynomials. A recent scheme for the numerical evaluation of Laguerre polynomials of any degree is described in [11].

As an example, the type of expansions in this paper have the following form. For the monic Laguerre polynomial of degree , we obtain:


This is expression (4.1) of the paper, specified to the standard associated Laguerre weight . It is valid for , where is related to through the MRS number , i.e., . The expansion itself follows from substituting the expansion of the matrix function . The leading order term is obtained from the identity matrix , and further terms are listed explicitly in Appendix A. We provide formulas and their implementation for an arbitrary number of terms in the asymptotic expansions, and compute up to terms in seconds on the architecture mentioned in § 7.4.

These formulas can also be applied to obtain asymptotic expansions of orthogonal polynomials with Hermite-type weights of the form on . In § 7.2, we show that they can be given in terms of asymptotics of Laguerre-type polynomials with , evaluated in .

We aim for a general non-polynomial weight function , though our results in this case are thus far not rigorously valid. In particular, we do not provide estimates for the remainder term. We do provide numerical indications that the expansions converge at the expected rate for increasing . Inspired by the requirements for the Jacobi case [6] and the technical conditions on in [17, §1], we conjecture that the expansions in this paper are valid as long as is analytic within the contours defined further on and grows faster than powers of for .

The structure of the paper is as follows. In § 2, we connect the Riemann-Hilbert problem for orthogonal polynomials that is analyzed in [28] with the expansion of a matrix-valued function and introduce some notation. We detail the Mhaskar-Rakhmanov-Saff (MRS) numbers and their asymptotic expansions for large in §3. The formulas for the asymptotic expansions of the polynomials in the different regions of the complex plane are stated in §4. We explain the computation of higher order terms of in §5 and provide a non-recursive definition for . Details on obtaining explicit expressions for higher order terms are provided in §6. We conclude the paper with a number of examples and numerical results in §7.

2 Asymptotic expansions for Laguerre–type polynomials

The largest root of a Laguerre-type polynomial grows with the degree . For example, it asymptotically behaves as for the standard associated Laguerre polynomials, with the (negative) zero of the Airy function closest to zero, see [19, 20, (18.16.14)] and [22, (6.32.4)]. The first step in the description of the asymptotics is to rescale the polynomials, such that the support of the zero-counting measure maps to the interval . The scaling is linear but -dependent and given by


where is the Mhaskar-Rakhmanov-Saff (MRS) number [17] defined further on in (3.1).

There is a distinction between several regions in the complex -plane, shown in Figure 1:

  • a complex neighbourhood of the interval excluding the endpoints, subsequently called the ‘lens’ (region \Romanbar1)

  • two disks around the endpoints and , called the right and left disk (regions \Romanbar3 and \Romanbar4)

  • and the remainder of the complex plane, the ‘outer region’ (region \Romanbar2).








Figure 1: Regions of the complex plane in which the polynomials have different asymptotic expansions, after rescaling the support of the zero-counting measure to the interval : the lens (\Romanbar1, a complex neighbourhood of the interval excluding the endpoints), the outer region (\Romanbar2, the remainder of the complex plane) and the right and left disks (\Romanbar3 and \Romanbar4, two disks around the endpoints and ).

2.1 Riemann–Hilbert formulation and steepest descent analysis

In this section, we briefly summarize the main features of the derivation in [28]. The approach is based on the Riemann–Hilbert formulation for orthogonal polynomials [10]: we seek a complex matrix-valued function that satisfies the following Riemann–Hilbert problem (RHP), cf. [28, §3]:

  1. : is analytic.

  2. has continuous boundary values , when going from the upper half-plane through the interval to the lower half-plane, respectively. These boundary values are related via a jump matrix:

  3. .

  4. The behaviour as is also specified, see [28, (3.3)].

It is proved in [10, 15], that the unique solution of this Riemann–Hilbert problem is

Here, the entry is the monic orthogonal polynomial. The entry relates to the polynomial of degree , while the second column contains the Cauchy transforms of both of these polynomials. Note that the weight function of the orthogonal polynomials enters through the jump condition in (b).

In order to obtain the large asymptotic behavior of , the Riemann–Hilbert formulation is combined with the Deift–Zhou steepest descent method for Riemann–Hilbert problems [8, 9]. In this case, the steepest descent analysis presented in [28] consists of the following sequence of (explicit and invertible) transformations:

These steps have a well-defined interpretation:

  • The first step is a normalization at infinity, such that


    This step comes at the cost of introducing rapidly oscillating entries in the new jump matrix for the Riemann-Hilbert problem for .

  • The second step is the opening of the so-called lens around : it factorizes the previous jump matrix such that outside of the lens in Figure 1, while is exponentially close to in in the upper and lower part of the lens. The shape of the lens is such that the oscillating entries on the diagonal in the jump matrix are transformed into exponentially decaying off-diagonal entries.

  • Finally, the last transformation gives rise to the disks in Figure 1 and is defined as


    Here, a global parametrix is introduced and constructed using the Szegő function (to be defined below). and are local parametrices that follow from a rather involved local analysis around the endpoints. We omit the details, but we note that is given explicitly in terms of standard Bessel and Hankel functions and their derivatives. The precise choice of these functions is made in such a way that satisfies a matching condition with the global parametrix on the boundary of the left disk, namely

    where will be defined in § 5.1 as the jump matrix for . A similar construction yields explicit expressions for in terms of the Airy function and its derivative. Since we know , and explicitly, we can determine the asymptotic expansion of in a closed formula.

The key idea is that the Riemann–Hilbert problem for can be solved explicitly in an asymptotic sense for large : it can be deduced that the matrix is itself close to the identity

uniformly for . Here, is a contour that results from the sequence of transformations outlined before, and consists of the boundaries of the regions in Figure 1. If we match all powers of and in creftype 2.3 via , we obtain higher order terms in the asymptotic expansions, which is exactly the technique outlined in § 5. Finally, reversing these transformations (since the different RHP are equivalent), one can obtain asymptotic information for as in different sectors of the complex plane, and in particular of the entry.

2.2 The function in the complex plane

The function is a matrix complex–valued function, analytic (element-wise) in , where consists of the boundaries of the regions in Figure 1. Also, for . Finally, as for polynomial with degree and independent of , there exist functions such that the function admits an asymptotic expansion of the form


We obtain different expressions for depending on the region in which lies. We will write and to refer to the coefficients for near and respectively, and to indicate the coefficients for outside these two disks. Our formulas for the asymptotic expansions are written in terms of these functions. One may simply substitute to obtain the leading order behaviour of the expansion. Higher-order expansions are obtained via recursive computation of the in §5, or alternatively using the explicit expressions listed in §A.

2.3 Auxiliary functions

We recall some terminology and notation from [28]. In the formulation of our results we use the third Pauli matrix and define for :

The following values arise in the recursive construction of the MRS numbers:

and we will also use the Pochhammer symbol or rising factorial

We will define the MRS number and the associated quantities and in § 3.2. The corresponding scaling (2.1) gives rise to the rescaled field :

We also define the following functions:


The function is non-standard in literature on asymptotics, but it is introduced here because it allows the statement of analytic continuations of some functions using standard branch cuts. We assume standard branch cuts of all analytic functions in this paper, such that the formulas are easily implemented. One may call , and phase functions for the orthonormal polynomials. They specify the oscillatory behaviour of for respectively away from the endpoints, near and near . Here, too, one can avoid specifying select branch cuts by not simplifying the definition of . The function corresponds to in [28] and is used for the analytic continuation of the polynomial in the left disk.

The conformal map from onto the exterior of the unit circle is used in the global parametrix , which determines the behaviour of away from and :


Finally, the coefficients and appear in asymptotics of Airy and modified Bessel functions in the local parametrices:

3 MRS numbers and related functions

The Mhaskar-Rakhmanov-Saff numbers satisfy [28]


We will explain how to compute these and quantities dependent on them for various types of : monomials, more general polynomials and more general analytic functions.

3.1 Monomial

If is monomial (), then [28]

We also recall from [28] the coefficients and polynomials with a slight adjustment for :


In the classical Laguerre case where , we have

The latter value for is well-known and it implies that the largest root of the Laguerre polynomial of degree grows approximately like .

3.2 General polynomial

For general polynomial , has an asymptotic expansion with fractional powers,


To compute the coefficients , we start from the equation at the end of the proof of [28, Prop 3.4]:


where we have defined


One uses the result [28, (3.8)]


and then recursively computes (3.5) for . Next, creftype 3.4 leads to

for and so on. We see that does not influence the MRS number, since it only rescales the weight function. The construction of from this section satisfies the condition creftype 3.1 asymptotically up to the correct order and also the following explicit result from [28, (3.8)]:


With these results, we can compute the polynomials and the coefficients as [28, §3]:


3.3 General function

For the calculation of the functions in case is not a polynomial, we introduce a numerical method. We provide an initial guess that satisfies to an iterative numerical procedure, so . The procedure finds a that approximately satisfies creftype 3.1, where the integral is computed by numerical integration. If needed, and can be approximated numerically as well. The other functions are given by integrals, [28, (3.11-16-38-40-24)]

The contour for should enclose the interval and the point . We choose to be a circle with a center halfway between the interval and , while still including and the interval. The integrals for and are also calculated numerically, so the former is computed by a double numerical integral.

Remark 3.1.

These expressions are also valid for polynomial . However, following the reasoning in this subsection only leads to a numerical value of for a given , as opposed to a full asymptotic expansion of in fractional powers of . The same observation holds for the functions defined above. In this case, the powers are implicitly present in all quantities that involve , while the results for polynomial are more explicit. We compare both approaches further in § 7.4 and 6.1.

3.4 Explicit expressions satisfied by the MRS numbers

Thus far we have obtained either asymptotic expansions of or a numerical estimation. The cases in which explicit expressions can be derived are limited, but in this section we aim to provide some more helpful expressions.

Inspired by creftype 2.7, we invert that conformal map by changing the coordinates in integral creftype 3.1:


The contour is half the unit circle, starting at through to . Note that is real-valued for on this halfcircle in the upper half of the complex plane. Hence, it is also real-valued when we take the complex conjugate of , corresponding to on the halfcircle in the lower half of the complex plane. If we assume that is real for real arguments, then the integral on the negative halfcircle, i.e. from through to , is the complex conjugate of the integral above. Combining both, we find that

where is a circle enclosing the origin in the counterclockwise direction.

The described change of variables maps a point in the interior of the unit circle to . If we is not entire, we need to substract additional residues; else, the contour encloses a single pole at the origin. In the case where is a polynomial of degree , we obtain from the residue theorem that is the root of a polynomial of degree :


We can remark that creftype 3.3 gives the asymptotic expansion of the zero of the -th degree polynomial creftype 3.12 in with respect to a factor in its constant coefficient. Exact solutions for are only available up to . For , this boils down to the standard associated Laguerre case . For we take the positive solution which also corresponds to creftypeplural 3.7 and 3.6:


We do find an explicit result for the non-polynomial function . In that case, we have


where denotes the Lambert-W function, and and are modified Bessel functions. For general , a similar technique may allow one to find an explicit expression satisfied by like creftype 3.15 without integrals. Solving that expression numerically avoids having to evaluate the integral creftype 3.1. However, it might become quite involved to derive higher order terms as explicitly as in § 6.4 and 6.3 from the resulting expansion of as .

4 Asymptotics of orthonormal polynomials and related coefficients

4.1 Lens \Romanbar1

Putting together the consecutive transformations in [28], for \Romanbar1 in Figure 1 and and related as in , we obtain


The asymptotics of are given in § 4.5. The full asymptotic expansion of is obtained by substituting the expansion for that we derive later on.

The asymptotic expansions of the orthonormal polynomials all separate two oscillatory terms (phase functions multiplied by ) from the non-oscillatory higher order terms. For polynomial of degree , the asymptotic expansion truncated after terms correspond to a relative error of size . In the special case of a monomial , the relative error improves to .

4.2 Outer region \Romanbar2

For \Romanbar2, the asymptotic expansion is


It may appear to be problematic that appears in the asymptotic expansions of the polynomials in the complex plane, especially for this region, since this factor grows very quickly. However, one may verify that this exponential behaviour is canceled out with other terms. More specifically, the term in (4.2) ensures that , .

4.3 Right disk \Romanbar3

The polynomials behave like an Airy function near the right endpoint (). This is typical asymptotic behaviour near a so-called ‘soft edge’, in the language of random matrix theory. Note that the in the following expression removes the branch cut, so that it can be used throughout , away from :


We would like to note a possible issue when computing the zeros of this expression, as one would do for Gaussian quadrature. The largest root for the standard associated Laguerre polynomials asymptotically behaves as , with the (negative) zero of the Airy function closest to zero, see [19, 20, (18.16.14)] and [22, (6.32.4)]. For a fixed but relatively high , this point may lie outside the support of the equilibrium measure [28, Rem. 3.8]. There will always be a larger for which the point lies inside. Still, for large one may want to pursue a different kind of asymptotic expansion, for example using asymptotics with a varying weight , as studied in for example [3], and apply it to the fixed .

4.4 Left disk \Romanbar4

The polynomials behave like a Bessel function of order near the left endpoint . For \Romanbar4, we obtain


It is not immediately obvious that the expansions (4.3) and (4.4) are analytic in the points and respectively. This will follow from the expression creftype 5.8 for and by also making a series expansion of the other terms at those points. For numerical purposes, it may be better to use those series expansions when evaluating close to (or at) and .

4.5 Asymptotics of leading order coefficients

The leading order coefficient of the orthonormal polynomials is , i.e. we have

where is the monic orthogonal polynomial of degree . For a monomial or more general function , the asymptotic expansion of is


The quantities are defined and extensively described in §5. They are the constant matrices that multiply and in the expansion for , of which we use the lower left elements here. Explicit expressions for these matrices up to are given in Appendix A. The constant coefficient only changes the scaling of the weight function and does not influence nor the matrices. However, it does influence through the coefficient , giving .

For general polynomial , the power of changes from to :