# Symmetric matrix formula for Gauss-Legendre quadrature

While searching the web, I came across the following algorithm for the Gauss-Legendre quadrature. I wasn't able to find a reference or a proof of my own as to why it works. I'll present it, and the question is - Why does it work?

The Method:

1. We construct the Matrix $A \in M_{N\times N} (\mathbb{R} )$ by $\forall 1\leq j\leq N-1 \, A_{j,j+1} = A_{j+1,j} = \frac{j}{\sqrt{4j^2 -1}}$, and all other entries are zero.

2. For some reason, the Legendre polynomial of order $N$ is the characteristic polynomial $p_N(x) = {\rm det}(x{\rm Id - A})$. Therefore, the eigenvalues of $A$ are the abscissas of the quadrature.

3. For the weights, we take $v^j$ to be the eigenvector of the $j-th$ smallest (negative) eigenvalue, and we have that $w_j =2(v^j _1) ^2, \, \forall 1\leq j \leq N$, with $v^j_1$ denoting the first coordinate.

I don't understand this method, nor can I find a proper reference to it in standard Numerical Integration books (Davis and Rebinowitz) or Orthogonal Polynomials books (Szego).

Meta - Note: in my old question I asked for a stable and efficient algorithm for the Gauß-Legendre quadrature, but none of the proposed methods was the one I've found. This is why I've opened a new question.

This is a particular implementation of a more general method, described in John Boyd's Why Eigenvalues Are Roots: A Derivation of the One-Dimensional Companion Matrix for General Orthogonal Polynomials (restricted access). Gauss-Legendre quadrature of order $n$ needs the roots of the Legendre polynomial $P_n(x)$ and a numerical root solver must guarantee that the roots are real. By reformulating the root finding problem into an eigenvalue problem for a symmetric matrix a numerical instability for large $n$ is avoided. The specific choice of symmetric companion matrix used here is derived on page 9 of this paper.

For the weights, I think the method described in the OP needs correction: It is not the $j$-th component of the eigenvector of the smallest eigenvalue that determines the weight $w_j$, but the first component of the $j$-th eigenvector. (I looked at the code linked by the OP, and it seems that is indeed what it does.) The relation is explained on page 223 of the Golub-Welsch paper, where all of this originates (or in more detail on page 6 of these notes).

• Thanks, you are right, I've corrected the formula. I'll sit down and read the references. – Amir Sagiv Apr 13 '16 at 13:16
• I can't seem to access Boyd's book or have it in our library. Is it available anywhere else? Is the same material covered somewhere else? – Amir Sagiv Apr 13 '16 at 13:32
• added a reference that derives the specific companion matrix used in the code, I think that paper is open access. – Carlo Beenakker Apr 13 '16 at 14:11
• Yep, this is exactly what I was looking for. Thanks! – Amir Sagiv Apr 13 '16 at 18:36

You are describing the Golub-Welsch algorithm. As described, it costs O(n^2) operations to compute n-point Gauss-Legendre quadrature rule.
http://www.ams.org/journals/mcom/1969-23-106/S0025-5718-69-99647-1/

If you are interested in large Gauss-Legendre quadrature rule, then better (faster and more accurate) algorithms now exist. See your previous question: Computing Gauss Legendre Quadrature for Large N