Computing Gaussian quadrature rules with high relative accuracy

The computation of n-point Gaussian quadrature rules for symmetric weight functions is considered in this paper. It is shown that the nodes and the weights of the Gaussian quadrature rule can be retrieved from the singular value decomposition of a bidiagonal matrix of size n/2. The proposed numerical method allows to compute the nodes with high relative accuracy and a computational complexity of O(n2).\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$ \mathcal {O} (n^{2}). $\end{document} We also describe an algorithm for computing the weights of a generic Gaussian quadrature rule with high relative accuracy. Numerical examples show the effectiveness of the proposed approach.


Introduction
The Golub and Welsch algorithm [12] is the classical way to compute the nodes λ i and the weights ω i of an n-point Gaussian quadrature rule (GQR), commonly used to approximate integrals of the following kind, where ω(x) ≥ 0 is a positive weight function in the interval [−τ, τ ], with τ > 0 and f (x) is a continuous function in the same interval. It computes the nodes, also called knots, as the eigenvalues of a symmetric tridiagonal matrix of order n, called the Jacobi matrix, while the corresponding weights are easily obtained from the first components of the associated eigenvectors (see [12] for more details). The nonzero entries of the Jacobi matrix are computed from the coefficients of the three-term recurrence relation of the orthogonal polynomials associated with the weight function ω(x).
The Matlab function gauss is a straightforward implementation of the Golub and Welsch algorithm [12], included in the package "OPQ: A Matlab suite of programs for generating orthogonal polynomials and related quadrature rules" by Walter Gautschi, and it is available at https://www.cs.purdue.edu/archives/2002/wxg/codes/ OPQ.html.
In [19] it is shown that, when the weight function ω(x) is symmetric with respect to the origin, the positive nodes of the n-point GQR are the eigenvalues of a symmetric tridiagonal matrix, of order n 2 if n is even and n+1 2 if n is odd, obtained by eliminating half of the unknowns from the original eigenvalue problem. Moreover, the Cholesky factor of the latter tridiagonal matrix of reduced order is a nonnegative lower bidiagonal matrix of dimensions n 2 × n 2 , if n is even, and n+1 2 × n−1 2 , if n is odd, whose entries in its main diagonal and lower subdiagonal are respectively the odd and the even elements of the first subdiagonal of the original Jacobi matrix.
In the same paper, the authors propose modifications of the Golub and Welsch algorithm for computing the nodes and weights of a symmetric GQR (SGQR), by means of a three-term recurrence relation corresponding to the orthogonal polynomials associated with the modified weight function. This resulted in faster Matlab functions than gauss.
Inspired by [19], in this paper we propose to compute the positive nodes of the npoint SGQR as the singular values of a nonnegative bidiagonal matrix. These can be computed in O(n 2 ) floating point operations and with high relative accuracy by the algorithm described in [5,6]. Moreover, the weights of the SGQR can be obtained from the first row of the right singular vector matrix. The stability of the three-term recurrence relation arising in computing the weights of a nonsymmetric n-point GQR is also analyzed and a novel numerical method for computing the weights with high relative accuracy is proposed.
A different approach for computing nodes and weights of a GQR in O(n) operations has been introduced in [10] and is implemented in the Matlab package Chebfun [15]. This method is based on approximations of the nodes and weights that appear to be relatively accurate in most nodes but no proof of such a result is given. For a very large number of nodes, this method is today the most efficient one [14,25].
The paper is organized as follows. The basic notation and definitions used in the paper are listed in Section 2. In Section 3, the main features of the Golub and Welsch algorithm are described. In Section 4, the proposed algorithm for computing the nodes of an n-point SGQR as the singular values of a bidiagonal matrix is described. Different techniques for computing the weights of an n-point GQR are described in Section 5. In Section 6 we show the elementwise relative accuracy of the weights computed by our new method and in Section 7, we give a number of numerical examples confirming our stability results. We then end with a section of concluding remarks.

Notation and definitions
Matrices are denoted by upper-case letters X, Y, . . . , , , . . . ; vectors with bold lower-case letters x, y, . . . , λ, θ, . . . ; scalars with lower-case letters x, y, . . . , λ, θ, . . .. The element i, j of a matrix A is denoted by a ij and the ith element of a vector x is denoted by x i .
Submatrices are denoted by the colon notation of Matlab: A(i : j, k : l) denotes the submatrix of A formed by the intersection of rows i to j and columns k to l, and A(i : j, :) and A(:, k : l) denote the rows of A from i to j and the columns of A and from k to l, respectively. Sometimes, A(i : j, k : l) is also denoted by A i:j,k:l and x(i : j) by x i:j .
The identity matrix of order n is denoted by I n , and its ith column, i = 1, . . . , n, i.e., the ith vector of the canonical basis of R n , is denoted by e i . The notation y stands for the largest integer not exceeding y ∈ R + . Given a vector x ∈ R n , then diag(x) is a diagonal matrix of order n, with the elements of vector x on the main diagonal. Given a matrix A ∈ R m×n , then diag(A) is a column vector of length i = min(m, n), whose entries are those of the main diagonal of A.

Computing the zeros of orthogonal polynomials
= 0, 1, . . . , be the sequence of orthonormal polynomials with respect to a positive weight function ω(x) in the interval where k is chosen to be positive. The polynomials p (x) satisfy the following threeterm recurrence relation [23] where Using (1), we can write the n-step recurrence relation The matrix J is called the Jacobi matrix [7]. The following theorem was shown in [7, Th. 1.31].
Theorem 1 Let J = Q Q T be the spectral decomposition of J, where ∈ R n×n is a diagonal matrix and λ := diag( ) = [λ 1 , . . . , λ n ] T , λ 1 < λ 2 < · · · < λ n , and Q ∈ R n×n is an orthogonal matrix, i.e., Q T Q = I n . Then p n (λ i ) = 0, i = 1, . . . , n, and Q = V D, with The eigenvalues λ i are the nodes of the GQR and the corresponding weights ω i are defined as (see [23]) Hence, the weights ω i of the GQR can be determined by (4), computing p (λ i ), = 0, 1, . . . , n − 1, i = 1, . . . , n, either by means of the three-term recurrence relation (1), as done in [19], or by computing the whole eigenvector matrix Q of J in (3). Both approaches will be described in Section 5, and their stability analysis will be provided as well. As shown in [26], the weights can be also obtained by the first row of Q as The Golub and Welsch algorithm [12], relying on a modification of the QR-method [4], yields the nodes and the weights of the GQR by computing only the eigenvalues of the Jacobi matrix and the first row of the associated eigenvector matrix.
Therefore, for weight functions symmetric with respect to the origin, it is sufficient to compute only the positive nodes and the corresponding weights. Without loss of generality, in the sequel we will consider the following reordering of λ, and V , i.e., Hence, the first n 2 entries of λ are the strictly positive eigenvalues of J in a decreasing order, the second n 2 entries of λ are the negative eigenvalues of J in an increasing order, and, if n is odd, the last entry of λ is the zero eigenvalue.

Computation of the nodes and the weights of a symmetric Gaussian quadrature rule
In this section we propose an alternative method to compute the positive nodes and corresponding weights of a SGQR. We will need the following Lemma [11].
Lemma 1 Let A ∈ R m×n , m ≤ n, and The following Corollary holds.

Corollary 1 Let B n be the bidiagonal matrix
with singular value decomposition where U and V are orthogonal, σ = σ 1 , σ 2 , . . . , σ n 2 T , and andˆ is the square diagonal matrix where P := I n (j , :) is the even-odd permutation matrix corresponding to the evenodd permutation of the index set i = [1, . . . , n], i.e., Proof Obviously, we have and (8) follows then straightforwardly from Lemma 1.
Taking the even-odd structure of P into account, we have e T 1 P T = e T n 2 +1 and thus the first row of Q is given by and, hence, by (4) and (5), the weights ω i are given by and, for n odd, Remark 1 Transforming the problem of computing the eigenvalues of J into the problem of computing the singular values of the nonnegative bidiagonal matrix B, allows us to compute the nodes of a SGQR with high relative accuracy by using the algorithm described in [5,6] and implemented in LAPACK (subroutine dlasq1.f) [1]. Note that dlasq1.f computes the singular values of square bidiagonal matrices and for the n-point SGQR with odd n, the bidiagonal matrix B is rectangular of size n 2 × n+1 2 . But one step of the QR iteration with zero shift, implemented as described in [5], is sufficient to transform B to B 0 , withB square bidiagonal of order n 2 to high relative accuracy. On the other hand, the componentwise stability of computing ω i , i = 1, . . . , n, from the first row of either Q or V , is not ensured if these matrices are computed by the QR algorithm. Indeed, these orthogonal matrices are computed in a normwise backward stable manner but not necessarily in a componentwise stable manner [27, p. 236], [16]. In Section 7, we show in the example of Fig. 2 that using the eigenvalue decomposition of J for computing the sequence (10), the computation of the small weights is not elementwise stable even though the eigenvectors are computed in a normwise backward stable manner [16]. For the sake of completeness, it is worth mentioning that a similar behavior is not observed for other classes of orthogonal polynomials, such as the Jacobi ones, since the weights do not have such a large range of values and, thus, the computation of GQR for such weights with the function gauss is as accurate as the method proposed in this paper.

Computation of the weights of the Gaussian quadrature rules
In this section we consider different techniques for computing the weights ω i , i = 1, . . . , n, of an n-point GQR, relying only on the knowledge of the corresponding nodes λ i for nonsymmetric weight functions. At the end of the section we will shortly describe how these techniques can be adapted to symmetric weight functions. As pointed out in Section 3, the nodes and the weights of the n-point Gaussian rule associated with a weight function ω(x), are the zeros λ i , i = 1, . . . , n, of the orthogonal polynomial p n (x) of degree n and the corresponding weights are The following methods have been proposed in the literature for computing the sequence 1. the eigenvalue decomposition [11]; 2. the forward three-term recurrence relation (FTTR) [8]; 3. the backward three-term recurrence relation (BTTR) [20,24].
The computation of the weights by means of the FTTR was considered in [19] without providing any stability analysis.
The stability of the generic FTTR and BTTR, i.e., sequences generated by threeterm recurrence relations not linked to orthogonal polynomials, was analyzed in [8,20,24], and in these papers the preference went to the BTTR for stability purposes.
For the weight functions depending on the parameters α and β listed in Table 1, we carried out extensive tests, choosing many different values of α and β, for computing the sequences p j (λ i ), j = 0, 1, . . . , n − 1, by means of FTTR and BTTR. FTTR works in an accurate way for some classes of weights and BTTR works in an accurate way for other classes of weights.
The results of these experiments are displayed in Table 1 in which we mention "Ok" if the method computes the sequence in an accurate way for all the considered values of α and β, and " " if the method fails to compute the sequence in an accurate way for some values of α and β. In Fig. 1 (left) we show on a logarithmic scale, the results obtained for the Hermite polynomials of degree j , j = 0, 1, . . . , 127, in the largest zero of the Hermite polynomial of degree 128, with FTTR in extended precision (128 digits), with FTTR in double precision and with BTTR in double precision. While the Hermite polynomials are accurately computed with FTTR, the results obtained by BTTR after some steps diverge from the actual values. In this case we Table 1 Success ("Ok") and failure ("No") results obtained by FTTR and BTTR when computing the sequence of orthogonal polynomials associated to the weight ω, considering many different values of α and β report an "Ok" in column FTTR and a " " in column BTTR in Table 1. In Fig. 1 (right) we show the results for the Hahn polynomials of degree j = 0, 1, 2, ..., 127, evaluated in the fourth smallest zero of the Hahn polynomial of degree 128, computed in double precision with FTTR and BTTR, for α = β = −.5. In this case, FTTR does not compute the sequence of Hahn polynomials in an accurate way, while BTTR does. We thus report a " " in column FTTR and an "Ok" in column BTTR in Table 1.
The Hahn polynomials in Table 1 are discrete orthogonal polynomials on the discrete set of n points {0, 1, 2, ..., n − 1} with respect to the weight α+k k β+n−1−k n−1−k , α, β > −1, k = 0, 1, 2, ..., n − 1. Hence, since experimentally it is not clear which of FTTR and BTTR can be chosen, we describe an algorithm, called LMV, that combines the FTTR and the BTTR methods in order to compute the sequence (10) with high relative accuracy. Letλ be a zero of p n (x) and let us denote byp and byp the following vectors, computed by FTTR and BTTR, respectively: . . .
The following theorem emphasizes the relationship between the sequencep j (λ), j = 1, . . . , n − 1, computed by (1) and the columns of theQ factor of the QR factorization of J −λI n .
since J −λI n is a tridiagonal matrix. Therefore, by Theorem 2, the subvector made by the first i entries of p(λ) is parallel to the vectors made by the first i entries of the columns j ofQ, with j ≥ i.
Theorem 3 shows the relationship between the sequence {p j (λ)} n−1 j =0 and the columns of the factorQ of the QL factorization of J −λI n . For the sake of brevity, we omit the proof since it is very similar to the one of Theorem 2.
withν i ∈ R,ν i = 0. Theorem 2 shows that the vectorp, computed using FTTR, can also be obtained by applying either one step of the implicit QR algorithm with shiftλ to J or computing the QR factorization of J −λI n , since the orthogonal matrices generated by both methods are the same. The last column of the orthogonal matrices generated by both methods will be parallel top. Therefore, forward instability can occur in the computation ofp if premature convergence occurs in one step of the forward implicit QR (FIQR) method with shiftλ to J [17,21].
On the other hand, Theorem 3 shows that the vectorp, computed applying BTTR, can also be obtained applying either one step of the backward implicit QR algorithm with shiftλ to J or computing the QL factorization of J −λI n . The first column of the orthogonal matrices generated by both methods will be parallel top. Therefore, forward instability can occur in the computation ofp if premature convergence occurs in one step of the backward implicit QL (BIQL) method with shiftλ to J [17,21].
The premature convergence of the implicit QR method with shiftλ depends on the distance between the eigenvalues of the consecutive matrices J 1:i,1:i and J 1:i+1,1:i+1 , i = 1, 2, . . . , n − 1. Hence, if λ (i) j , j = 1, . . . , i, and λ (i+1) j , j = 1, . . . , i + 1, are the eigenvalues of J 1:i,1:i and J 1:i+1,1:i+1 , respectively, then, by the Cauchy interlacing Theorem [22], j = 1, . . . , i, i = 1, . . . , n − 1, are sufficiently large, then premature convergence does not occur in one step of the implicit forward QR algorithm with shiftλ and, hence, FTTR computes the sequence (10) accurately. This is the reason why FTTR computes the sequence (10) accurately for the Chebyshev polynomials T j (x) of the first kind, i.e., Jacobi polynomials associated with the weight 1 √ 1−x 2 in the interval [−1, 1], since the distance of the zeros of two consecutive polynomials T j (x) an T j +1 (x) is of order O 1 j 2 at least. In [17,18], an algorithm combining one step of FIQR and one step of BIQL with shiftλ is described in order to compute the corresponding eigenvector. Therefore, the eigenvector associated withλ is computed combining the firstj − 1 Givens rotations of FIQR with shiftλ and the first n −j rotations of BIQL with shiftλ. It is proven that each eigenvector is computed accurately with O(n) floating point operations. Once the eigenvector is computed, the weights are obtained by applying (9).
Following [17,18], we now describe a recursive procedure to determine an interval in which the indexj lies. Let us consider the sequence of Givens rotationsG i ∈ R n×n , i = 1, . . . , n − 1, defined in (11). It turns out that . . .
This suggests the Algorithm 1, written in a Matlab-style, to compute the interval [ĵ,j ] in whichj lies.
Once the interval [ĵ,j ] is determined, the indexj is chosen as the index with the maximum element in absolute value in the subvectorpĵ :j [17,18].

Remark 3
Givenλ, a similar recursion holds for estimating if premature convergence occurs in a singular value of the bidiagonal matrix B.

Remark 4
We have described in Section 4 that the sequence p (λ i ), = 0, 1, . . . , n − 1, can be retrieved from the eigenvector of J associated with λ i . Since the eigenvalue decomposition of J can be obtained from the singular value decomposition of B, then, the eigenvector sequences can be retrieved in a similar way from the left and right singular vectors associated with the singular value λ i of the corresponding bidiagonal matrix. For the sake of brevity, we omit the details.

Stability of the eigenvectors and weights
In this section we analyze the sensitivity of the calculation of the eigenvector of the tridiagonal matrix J and the sensitivity of the corresponding weight of the GQR, for a particular eigenvalue λ j . Let us define the shifted matrix as T := J − λ j I n , where T is tridiagonal and unreduced, let p be the corresponding vector of orthogonal polynomials evaluated at λ j , i.e., T p = 0, and let ω be the corresponding weight, i.e., ω := 1/(p T p). Let us now consider any nonzero element p i = 0 of p. We then partition the rows of the matrix T as follows: This then yields the following theorem.
and for any | q i |= q ∞ , or, equivalently, | p i |= p ∞ , we have Proof The Cauchy inequalities for singular values yields σ n−1 (T (i) ) ≤ σ n−1 (T ) since we deleted one row of the matrix T to obtain T (i) . Let Q ∈ R n×(n−1) be the orthogonal complement of q, i.e., Q T Q = I n−1 and Q T q = 0. Then Let v ∈ R n−1 be such that v 2 = 1, and MT Qv 2 = σ min (MT Q) = σ min (T (i) Q).
We then have MT Qv 2 ≥ σ min (M) T Qv 2 ≥ σ min (M)σ min (T Q), which implies Moreover, MM T = I + e iq T +qe T i , and it has been shown in [2,Thm 3.6] that Putting this together yields (16). The inequalities (17) for | q i |= q ∞ then follow from | q i |≥ 1/ √ n.
Let us denote by p f l the approximation of the vector p computed by LMV, and let us look at how p f l is constructed from the matrix equation T (i) p = 0 and show that it is backward stable in the sense that there exists a perturbation such that the computed vector p f l satisfies exactly the equation where M is the machine precision of the computer used. Once the index i has been chosen, the LMV method computes the vector shared by the kernels of T 1:i−1,: and T i+1:n . Basis vectors for these kernels are respectively given by wherex andx are arbitrary. In order to construct a common vector in the two kernels, we impose β = αp i /p i , and then choose α such that the common vector corresponds to the initializationp 0 (λ j ) = 1/ √ μ 0 . The subvectors α p 1:i−1 p i and β p î p i+1:n , are computed by a forward and backward recurrence using the three-term recurrence of the tridiagonal matrix T . After fixing the starting values in these recurrences, they each can be interpreted as a back substitution of a triangular system of equations. This was analyzed in depth in [16,Ch. 8], from which it follows that the computed We point out that the stability result for the eigenvector is similar to what one has for the eigendecomposition of the matrix J , since it is inversely proportional to the smallest nonzero singular value of T , i.e., to the smallest gap between λ j and the remaining eigenvalues. The sensitivity of each weight has this same inverse factor, but it has, except for this factor, a forward error that is stable in a relative sense. This is a strong property that is not shared by the eigenvalue method.

Numerical examples
In this section we compare the computation of the nodes and weights of n-point GQRs obtained by the proposed method, called LMV, to gauss, the Matlab function available in [9] and to the methods proposed in [19]. All the experiments were performed in Matlab ver. R2020b. In [19], the authors show that the positive nodes and the corresponding weights of a GQR associated with a symmetric weight function can be retrieved from the eigendecomposition of a tridiagonal matrix J n 2 of order n 2 , providing the Matlab functions displayed in Table 2. In the first and second example these methods are used to compute n-point GQRs corresponding to the Chebyshev weights of first and second kind, respectively, since their nodes and weights are known. In the third example the proposed method is compared to the function gauss in computing an integral on the whole real line.

Example 1
The nodes of the n-point GQR associated with the Chebyshev weight of the first kind, , are x j = cos( (2j −1)π 2n ), j = 1, . . . , n, the zeros of T n (x), the Chebyshev polynomial of the first kind of degree n. Moreover, the weights are w j = π/n, j = 1, . . . , n. The maxima of the relative errors of the nodes computed by the considered numerical methods  Table 3, while the maxima of the relative errors of the computed nodes Table 4. In all the cases, the relative errors of the nodes and weights computed by the proposed method are comparable to those of the results yielded by the algorithms proposed in [19].

Example 2
The nodes of the n-point GQR associated with the Chebyshev weight of the second kind, are x j = cos( π n+1 ), j = 1, . . . , n, the zeros of U n (x), the Chebyshev polynomial of the second kind of degree n. Moreover, the weights are w j = (1 − x 2 j ) π n+1 , j = 1, . . . , n. The maxima of the relative errors of the nodes computed by the considered numerical methods are reported in Table 5 while the maxima of the relative errors of the computed nodes are reported in Table 6.
In all the cases, the relative errors of the nodes and weights computed by the proposed method are smaller than those of the results yielded by the algorithms proposed in [19].

Example 3
In this example we consider a GQR for integrals on the whole real line with the Hermite weight ω(  We computed the Hermite weights, for different values of n, with the following three different methods: • the function gauss [9], computed in double precision. • the function gauss [9], computed in variable precision with 128 digits; these values can therefore be considered as exact values for the weights. • the proposed method, called LMV.    , computed with the new method LMV (in " "). The corresponding relative error estimate e est (ω LMV j ) for each weight computed with the method LMV is also given in "•" In this paper it is shown that, for symmetric weight functions, the positive nodes and the corresponding weights can be computed from a bidiagonal matrix with positive entries of size n 2 × n+1 2 . Therefore the nodes can be computed with high relative accuracy by an algorithm proposed by Demmel and Kahan. Moreover, the stability of different methods for computing the weights is analyzed, proposing an algorithm for computing them with relative accuracy of the order of the machine precision. The numerical experiments confirm the effectiveness of the proposed approach.