Generalized eigenvalue methods for Gaussian quadrature rules

A quadrature rule of a measure $\mu$ on the real line represents a convex combination of finitely many evaluations at points, called nodes, that agrees with integration against $\mu$ for all polynomials up to some fixed degree. In this paper, we present a bivariate polynomial whose roots parametrize the nodes of minimal quadrature rules for measures on the real line. We give two symmetric determinantal formulas for this polynomial, which translate the problem of finding the nodes to solving a generalized eigenvalue problem.


INTRODUCTION
Given a (positive Borel) measure µ on R, a classical problem in numerical analysis is to approximate the integral with respect to µ of a suitably well-behaved function f . One approach is via so called quadrature rules. These approximate the integral by a weighted sum of function values at specified points. One classical construction for quadrature rules designed to approximate the integral of continuous functions consists of demanding an exact evaluation of the integral for all polynomials of degree ≤ D. If the moments of µ exist and are finite, then this amounts to finding a measure supported on finitely-many points whose moments agree with those of µ up to degree D.
We use t as a formal variable on the real line and and write R[t] ≤D for the vector space of real polynomials of degree at most D. For k ∈ N 0 , we denote the kth moment of µ, if it exists and is finite, by m k = t k dµ.
Definition 1. Suppose D is a positive integer and µ is a measure on R whose moments up to degree D exist and are finite. For x ∈ R ∪ {∞}, define the linear function ev x : We sometimes write ev D ∞ for ev ∞ to emphasize its dependence on D. A quadrature rule of degree D for µ is a finite set N ⊂ R ∪ {∞} together with a function w : N → R >0 with We call N the nodes of the quadrature rule.

Remark 2.
In Definition 1 we allow nodes at infinity. This is not desirable for application in numerical analysis, since one cannot generally evaluate functions at these nodes. However, compactifying the real line by adding ∞ makes certain arguments easier. Although our statements are phrased with this more general notion of quadrature, our main result presented below actually provides a tool to explicitly distinguish the cases of quadrature rules with all nodes real from the case where one of the nodes is ∞. Furthermore, it is a classical theorem that quadrature rules for µ with no nodes at infinity exist for every degree D provided the moments of µ exist and are finite up to the same degree (see e.g. [La1,Theorem 5.8], [La2,Theorem 5.9], [Sch,Theorem 1.24]).
We say the measure µ is non-degenerate in degree d if its moments m k are finite up to degree k = 2d and for every nonzero, nonnegative polynomial f ∈ R[t] ≤2d we have f dµ > 0. Since in one variable any nonnegative polynomial f is sum of squares of polynomials, this is equivalent to demanding that p 2 dµ > 0, for every 0 = p ∈ R[t] ≤d . This property can be checked quite conveniently in the following way.
Definition 3. Consider the quadratic form p ∈ R[t] → p 2 dµ and restrict it to R[t] ≤d . With respect to the monomial basis 1, t, . . . , t d for R[t] ≤d , this quadratic form is represented by the (d + 1) × (d + 1) Hankel matrix with (i, j)th entry m i+j−2 : With this notation, a measure possessing finite moments up to degree 2d is non-degenerate in degree d if and only if det(M d ) = 0.
From the point of view of numerical analysis it is desirable to have a quadrature rule that is exact up to a certain degree and requires the fewest number of evaluations possible. Such a quadrature rule is called a Gaussian quadrature rule for µ.
It is known that when D = 2d + 1 is odd and µ is non-degenerate in degree d, there is a unique quadrature rule for µ with d + 1 nodes, and none with fewer nodes [Sch,Cor. 9.8].
The nodes can be found as follows. Let M d denote the (d + 1) × (d + 1) matrix representing the quadratic form p → t · p 2 dµ with respect to the monomial basis 1, t, . . . , t d of R[t] ≤d , meaning that (M d ) i,j equals m i+j−1 . Then the d + 1 nodes of the unique Gaussian quadrature rule for µ in degree 2d + 1 are the d + 1 roots of the univariate polynomial det(xM d − M d ), see e.g. [Sze,form. 2.2.9,p. 27]. Since M d and M d are real symmetric matrices and M d is positive definite, this writes the problem of finding these d + 1 nodes as a generalized eigenvalue problem.
In this paper we focus on the case when D = 2d is even. In this case, there is a oneparameter family of quadrature rules for µ with d + 1 nodes (see e.g. [Sch,Thm 9.7]). Here we reprove this fact by constructing a polynomial F ∈ R[x, y] with degree d in each of x and y and the property that for every y ∈ R, the d roots of F(x, y) ∈ R[x] ≤d are the other d nodes (among them possibly ∞) of this unique quadrature rule with y as a node. Furthermore, we give symmetric determinantal representations of F, which again translates the problem of finding nodes of a quadrature rule into finding the generalized eigenvalues of a real symmetric matrix, i.e. solving det(xA − B) = 0 where A, B are real symmetric matrices and A is positive semidefinite (see e.g. [BDDRV,Chapter 4]). While the literature on Gaussian quadrature rules is vast, to our knowledge these formulas are new.
To construct F and its determinantal representations, we consider three quadratic forms on R[t] ≤d−1 , namely those taking p ∈ R[t] ≤d−1 to p 2 dµ, t · p 2 dµ, and t 2 · p 2 dµ, respectively. We call M d−1 , M d−1 and M d−1 the matrices representing these quadratic forms with respect to the basis 1, t, . . . , t d−1 . That is, define the d × d matrices Theorem 4. Let µ be a measure on R that is non-degenerate in degree d ≥ 1. There is a unique (up to scaling) polynomial F ∈ R[x, y] of degree 2d with the property that for x, y ∈ R, F(x, y) = 0 if and only if there is a Gaussian quadrature rule for µ of degree 2d with nodes x = r 1 , y = r 2 , r 3 , . . . , r d+1 in R ∪ {∞}. This polynomial has the following two determinantal representations: (a) the determinant of the d × d matrix with bilinear entries in x, y, Moreover, for all x ∈ R, with det(xM d−1 − M d−1 ) = 0, all nodes are on the real line, i.e., the associated quadrature rule has no nodes at infinity.
We prove the two parts of this theorem in Sections 2 and 3. The implications for finding quadrature rules as generalized eigenvalues are made explicit in Remarks 6 and 13. In Section 4, we discuss possible generalizations of this result.

BILINEAR DETERMINANTAL REPRESENTATION
In this section, we prove Theorem 4(a). The proof relies heavily on understanding a particular symmetric bilinear form B x,y on R[t] ≤d−1 . For x, y ∈ R and p, q ∈ R[t] ≤d−1 , let Note that the symmetric matrix xyM d−1 − (x + y)M d−1 + M d−1 represents B x,y with respect to the basis 1, t, t 2 , . . . , t d−1 . Define F ∈ R[x, y] to be the polynomial We will show that F satisfies the requirements of Theorem 4.
Proof of Theorem 4(a). Suppose that there exists a quadrature rule for µ of degree 2d with nodes r 1 , r 2 , . . . , r d+1 in R ∪ {∞}, meaning that there exist w 1 , . . . , w d+1 ∈ R >0 so that Note that because µ is non-degenerate in degree d, the nodes r 1 , . . . , r d+1 are necessarily distinct, so that after reindexing we may assume that r 1 , . . . , r d ∈ R.
We claim that F(r 1 , r 2 ) = 0. To see this, let q be the unique (up to scaling) nonzero The last equality follows from (2) and the fact that deg(p) ≤ d − 1. Therefore q is an element of the right kernel of B r 1 ,r 2 . Since B r 1 ,r 2 drops rank, the determinant F(r 1 , r 2 ) of the representing matrix equals zero.
Conversely, suppose that for x, y ∈ R, F(x, y) = 0. Then the kernel of B x,y contains a polynomial q ∈ R First, we argue that x, y, and the roots of q are real and pairwise distinct and that the degree of q is d − 2 or d − 1. If not, there would exist a non-zero polynomial p ∈ R[t] ≤d−1 for which f = (t − x)(t − y)pq is nonnegative on R. Since f dµ = 0, this contradicts the assumption that µ is non-degenerate in degree d.
Let r 1 = x, r 2 = y and denote the roots of q by r 3 , . . . , r d+1 , where we take r d+1 = ∞ if deg(q) = d − 2. Consider the conic hull of the d + 1 points ev r 1 , . . . , ev r d+1 in R[t] * ≤2d . This is a (d + 1)-dimensional simplicial convex cone in the (2d + 1)-dimensional vectorspace R[t] * ≤2d . Therefore this cone is defined by d linear equalities and d + 1 linear inequalities, which we will now identify by inspection. For each i = 1, . . . , d + 1, let f i be the unique (up to scaling) nonzero polynomial of degree ≤ d for which ev r j ( f i ) = 0 for each j = i. ≤2d . It therefore follows that the conic hull of ev r 1 , . . . , ev r d+1 is the set of L ∈ R[t] * ≤2d satisfying Since µ is non-degenerate in degree d, each of these weights must be positive.
We remark that for every x ∈ R, the matrix M(x, x) is positive definite. To see this, note that M(x, x) represents the quadratic form on R[t] ≤d−1 given by p → p 2 (t − x) 2 dµ with respect to the monomial basis. Since µ is non-degenerate in degree d, this is positive for any 0 = p ∈ R[t] ≤d−1 .
Remark 6. For fixed y ∈ R, this allows us to find the roots of F(x, y) ∈ R[x] as the generalized eigenvalues of a d × d real symmetric matrix. If y is larger than all of the roots of det We can find the roots in x as the following generalized eigenvalue problem: is not positive definite, this formula still holds, but may not be as numerically stable. We can instead make a change of variables x =x + y, which gives In particular, this has the formxA + B where B is positive definite. Suppose λ is a root of det(λB − A). Thenx = −1/λ is a solution to det(xA + B) = 0. Therefore Note that this even works when λ = 0 if we take y − 1/λ to be ∞.

Corollary 7.
Let µ be a measure on R that is non-degenerate in degree d ≥ 1. For y ∈ R there is a unique quadrature rule for µ of degree 2d with d + 1 nodes, one of which is y.
Since the matrix pencil {M(x, y) : x ∈ R} contains a positive definite matrix M(y, y), the roots of F(x, y) = 0 must all be real. In particular the existence of one real root x implies, by Theorem 4, the existence of a quadrature rule of degree 2d of µ whose d + 1 nodes include y. Moreover the other d nodes x, r 3 , . . . , r d+1 are necessarily the d roots of F(x, y) = 0. As in the proof of Theorem 4, the conic hull of ev x , ev y , ev r 3 , . . . , ev r d+1 is a simplicial cone containing L µ , meaning that there is a unique representation of L µ as a nonnegative combination of them.
For fixed y ∈ R, the polynomial F(x, y) ∈ R[x] has three real roots, r 1 , r 2 , r 3 ∈ R ∪ {∞}, which, together with y, form the nodes of a quadrature rule for µ of degree 6. The matrix M(x, y) has the form xA + B for some real symmetric matrices A and B (depending on y). The roots of F(x, y) are the roots of det(xA + B).
Moreover, by making a change of coordinates we can make A positive definite. For example, for y = 1, we set x = 1/λ + 1 to get which has the form λA − B, where A = M(1, 1) is positive definite. Solving the generazlied eigenvalue problem det(λA − B) = 0 gives λ ≈ −0.66, −0.32, 0.60. The solutions to F(x, 1) = 0 are then x = 1 + 1/λ ≈ −2.15, −0.52, 2.67. Thus there is a quadrature rule for µ of degree 6 with nodes consisting of 1 and these three roots. The curve given by F(x, y) = 0 along with the line y = 1 are shown in Figure 1.
For y ∈ {0, ± √ 3}, the polynomial F(x, y) ∈ R[x] becomes quadratic with two real roots {0, ± √ 3}\{y}. From this we conclude that there is a quadrature rule of degree 6 for µ with the four nodes 0, ± √ 3, ∞. For y > √ 3, the matrix coefficient of x in M(x, y) is positive definite. In this case, M(x, y) already has the form xA − B where A is positive definite, so no change of coordinates is required to translate this into a generalized eigenvalue problem.

LINEAR DETERMINANTAL REPRESENTATION
In this section we prove Theorem 4(b). As in Section 2, we construct a bilinear form depending on x, y ∈ R and construct a non-zero kernel of it for those pairs x, y that can be extended to d + 1 nodes of a quadrature rule for µ of degree 2d.
For x, y ∈ R, we define a bilinear form B x,y on R ⊕ R[t] ≤d−1 ⊕ R[t] ≤d−1 ∼ = R 2d+1 . Given p = (p 0 , p 1 , p 2 ) and q = (q 0 , q 1 , Choosing the basis 1, t, . . . , t d−1 for both copies of R[t] ≤d−1 represents this symmetric bilinear form as the (2d + 1) × (2d + 1) matrix given in Theorem 4 where M d−1 and M d−1 are the d × d matrices defined in (1). In order to construct a non-zero element in the kernel of B x,y , we need to build up some basic facts. The first concerns quadrature rules for µ whose nodes include ∞. As this polynomial will be used heavily in the text below, denote Lemma 9. The polynomial F ∞ has d real roots s 1 , . . . , s d ∈ R and there exist weights w 1 , . . . , w d ∈ R >0 for which Proof. Since M d−1 is positive definite, F ∞ has d real roots s 1 , . . . , s d ∈ R, up to multiplicity. For any root s, the matrix sM d−1 − M d−1 has rank < d, meaning that the polynomial For any r ∈ R with F(r, s) = 0, there is a quadrature rule for µ of degree 2d with d + 1 nodes containing s, r, and ∞. Then the unique quadrature rule of degree 2d with d + 1 nodes containing r also contains the node ∞. This implies that F(x, r) has degree ≤ d − 1 and F ∞ (r) = 0, meaning that r ∈ {s 1 , . . . , s d }. Therefore there is a quadrature rule of degree 2d for µ with nodes s 1 , . . . , s d , ∞. In particular, there exist w 1 , . . . , w d , w ∞ ∈ R >0 for which (3) holds. Then w ∞ is the largest λ for which the quadratic form p → p 2 dµ − λ ev ∞ (p 2 ) is nonnegative on R[t] ≤d . This is the largest λ for which the matrix M d − λe d+1 e T d+1 is positive semidefinite. We find this by solving the equation det(M d − λe d+1 e T d+1 ) = 0, which gives λ = w ∞ = det(M d )/ det(M d−1 ). Lemma 10. Let w ∞ , s 1 , . . . , s d ∈ R be as given by Lemma 9. If y ∈ R with F ∞ (y) = 0, then there is a quadrature rule for µ of degree 2d with nodes y, r 1 , . . . , r d ∈ R. Let Then for all p ∈ R[t] ≤d−1 , Proof. By Corollary 7, there is a quadrature rule for µ of degree 2d with nodes y, r 1 , . . . , r d in R ∪ {∞}. Since F ∞ (y) = 0 the univariate polynomial F(x, y) ∈ R[x] has degree d and by the uniqueness of the quadrature rule, r 1 , . . . , r d are its roots. In particular, r j ∈ R. Then for p ∈ R[t] ≤d−1 , we have that . The first equality comes from the fact that there is a quadrature rule for µ of degree 2d with nodes y, r 1 , . . . , r d . The second follows from the equality in Lemma 9.
We now make a special choice of q and use Lemma 10 to greatly simplify B x,y (p, q).
Lemma 11. Suppose x, y ∈ R satisfy F ∞ (x) = 0 and F ∞ (y) = 0. Take w ∞ ∈ R and q x , q y ∈ R[t] ≤d−1 as given by Lemmas 9 and 10 and fix q = (w ∞ , q x , −q y ). Then for any p = (p 0 , p 1 , Proof. Let q = (w ∞ , q x , −q y ) and suppose p = (p 0 , p 1 , Using the properties of q x , q y given in Lemma 10, this simplifies to and have the same real roots along the line y = y for y ∈ R, contradicting the distinctness of these roots. Thus it suffices to show that a polynomial vanishing on the real variety of H must be a polynomial multiple of each H i . Since each H i satisfies our assumption as well, we can assume without loss of generality that H is irreducible itself. The real variety of H contains infinitely many points so its Zariski closure in C 2 is at least one dimensional. By irreducibility of H, its real variety Zariski-dense in its complex variety. Now the claim follows from Hilbert's Nullstellensatz. Now we are ready to prove Theorem 4(b).
Proof of Theorem 4(b). Let F be the polynomial given by Theorem 4(a). Then F has degree 2d in x and y, with top degree part det(M d−1 )x d y d , which is non-zero by the nondegenerateness of µ. Then (x − y)F has degree 2d + 1 with top degree part equal to det(M d−1 )(x − y)x d y d . Also, for every y ∈ R, the polynomial (x − y)F(x, y) ∈ R[x] has distinct, real roots. Lemma 12 then implies that any polynomial G ∈ R[x, y] vanishing on the real variety of (x − y)F must be a polynomial multiple of it.
We further claim that the points ( ) also vanishes on V R ((x − y)F) and is therefore a multiple of (x − y)F. It suffices to show that (x − y)F has no factors in common with F ∞ (x)F ∞ (y). The factors of F ∞ (x)F ∞ (y) are given by Lemma 9. Suppose for the sake of contradiction that (x − s i ) divides (x − y)F for some i = 1, . . . , d. It must be that (x − s i ) divides F. This implies that F(s i , s i ) equals zero, which contradicts the observation in Remark 5 that the matrix M(s i , s i ) is positive definite and F(s i , s i ) = det(M(s i , s i )) > 0.
Let G denote the determinant of the (2d + 1) × (2d + 1) matrix representing the bilinear form B x,y . We will show that ( Note that c · G is also a polynomial of degree 2d + 1. Inspection shows that its top degree part is det(M d−1 )(x − y)x d y d . Thus by the above argument, it suffices to show that G(x, y) = 0 for all (x, y) Now take x, y ∈ R with (x − y)F(x, y) = 0 and F ∞ (x)F ∞ (y) = 0. Let q x , q y ∈ R[t] ≤d−1 be the polynomials given by Lemma 10 and let q = (w ∞ , q x , −q y ). We claim that q is in the kernel of B x,y . To see this, let p = (p 0 , p 1 , by Lemma 11. If x = y, this is clearly zero. Otherwise F(x, y) = 0 and by Theorem 4(a) there exists r 3 , . . . , r d+1 ∈ R ∪ {∞} and a quadrature rule of degree 2d for µ with nodes x, y, r 3 , . . . , r d+1 . Since F ∞ (x) = 0, each r i ∈ R. Then Expanding and looking at the coefficient of t d−1 reveals that In particular, ev ∞ (q x − q y ) = y − x, giving B x,y (p, q) = 0. Since the bilinear form B x,y has a non-zero kernel, the determinant, G(x, y), of its representing matrix is zero.
Remark 13. Again this translates the problem of finding the roots of F(x, y) for fixed y ∈ R into a generalized eigenvalue problem. Consider the (2d + 1) × (2d + 1) symmetric matrix polynomial Since Then M(x, 1) has the form xA − B where A is a positive semidefinite matrix of rank four. The determinant det(M(x, 1)) has four roots, x ≈ −2.15, −0.52, 1, 2.67, which are the nodes of a quadrature rule for µ of degree 6. The curve det(M) = 0 along with the line y = 1 are shown in Figure 2.
Remark 15. Note that the matrix M(x, −y) has the form xA + yB + C where the matrices A, B are positive semidefinite. For any a 1 , a 2 ∈ R >0 and b 1 , b 2 ∈ R, the matrix coefficient of t in M(a 1 t + b 1 , −a 2 t + b 2 ) is positive definite, implying that the polynomial F(a 1 t + b 1 , −a 2 t + b 2 ) ∈ R[t] is real-rooted. One can see this in Figures 1 and 2, as any line with negative slope intersects the curve V(F) in six real points. It also shows the polynomial F(x, −y) is real stable [Wag,Prop. 2.4]. The Helton-Vinnikov theorem [HV,Thm. 2.2] then implies that not only (x − y) · F has a (2d + 1) × (2d + 1) linear determinantal representation like in the Theorem 4(b) but also F itself has a determinantal representation F = det(xA + yB + C), where A, B, C are 2d × 2d real symmetric matrices and A and −B are positive semidefinite. However it is unclear if there exists such a representation for which the entries of A, B, C are easily calculated from the moments m k of µ, as with M.

UNIVARIATE QUADRATURE RULES WITH MORE NODES
It is natural to try to generalize the above discussion to situations where more nodes of a quadrature rule are prescribed. Finding a quadrature rule means specifying an even number of real parameters since each node comes with a weight. We will now consider the following minimal problem: Problem 16. For integers n, ≥ 1, we are given n + 2 + 1 moments m 0 , . . . , m n+2 of a (positive Borel) measure µ on the real line that is non-degenerate in degree n + 2 and n − 1 real numbers x 1 , . . . , x n−1 . Does there exist a quadrature rule for µ of degree n + 2 with n + nodes including x 1 , . . . , x n−1 ?
Specifying the quadrature rule requires 2( + 1) + n − 1 = n + 2 + 1 parameters, as we have two parameters for each of + 1 unspecified nodes and one parameter for the weight of each of the n − 1 specified nodes. Therefore the number of parameters that we have to choose is exactly equal to the number of moments that we need to match.
For n = 1, the problem is solved by the well-known Gaussian quadrature. The case n = 2 is the main focus of this paper, and it is classically known that such a measure exists. Unfortunately, for n = 3 this can fail even if the measure is non-degenerate in every degree. To see this, we need some preparation.
For each integer 0 ≤ k ≤ n, consider the quadratic form p → t k · p 2 dµ on R[t] ≤ , which, with respect to the basis 1, t, . . . , t , is represented by the ( + 1) × ( + 1) matrix Notice that the highest moment of µ needed to specify these matrices is m n+2 , achieved by i = j = + 1 and k = n.
Therefore the bilinear form B X has a nonzero kernel and the determinant of its representing matrix (4) is zero.
Unfortunately, the converse of Proposition 17 does not hold for n > 2.
Example 18. Consider the measure given by the exponential distribution on R, whose k-th moment is m k = k! for every k ∈ N 0 . Let n = 3, = 3 so that n + 2 = 9. Then in Problem 16, we want to build a quadrature rule for µ of degree 9 with at most n + = 6 nodes including n − 1 = 2 specified nodes x 1 , x 2 ∈ R. However for x 1 = 1 3 and x 2 = 11, the determinant of the 4 × 4 matrix in Proposition 17 equals (up to a positive constant) 137503x 4 3 − 1695024x 3 3 + 11282760x 2 3 − 41197920x 3 + 46998216 ∈ R[x 3 ], which has complex roots x 3 ≈ 1.87, 5.20, 2.63 ± i 5.31, only two of which are real. Since µ is non-degenerate in every degree, any quadrature rule for µ in degree 9 has at ≥ 5 nodes. However, by Proposition 17 there are only two possibilities for the remaining ≥ 3 nodes. Therefore there is no quadrature rule for µ of degree 9 with ≤ 6 nodes including x 1 = 1 3 and x 2 = 11.