# Integrating Lagrange polynomials with many nodes, round-off

Given a set of points $\{x_j\}_{j=1}^n$ in $[-1, 1]$, I would like to compute $$\int_{-1}^{1} L_i(x)\,\text{d} x$$ exactly. $L_i$ is the Lagrange polynomial with respect to the points $x_j$ with $x_i$ as node, i.e., $$L_i(x) = \prod_{j\neq i} \frac{x - x_j}{x_i - x_j}.$$ Since this is a polynomial of degree $n$, I could use any old Gaussian quadrature of sufficient degree. This works well if $n$ isn't too large, but leads to results flawed by round-off errors for large $n$.

Any idea how to avoid those?

• This depends on where $x_j$'s are, but have you checked that your $L_i$'s are well-behaved? In the worst case, with $x_j$ being uniformly distributed, you get the Runge phenomenon ($L_i$'s oscillatory and large), in which case it's not really roundoff errors causing trouble. – Kirill Mar 16 '17 at 19:32
• Also, nitpick: dividing by small numbers is a well-conditioned operation, it is rather the subsequent subtraction of large nearly-equal numbers that is ill-conditioned and leads to numerical instability. – Kirill Mar 16 '17 at 19:33
• It seems like you're trying to compute $(2,0,\frac23,0,\frac25,0,\ldots)^\top V^{-1}$ where $V$ is the Vandermonde matrix of $x_j$'s. Can you say what the condition number of $V$ is? – Kirill Mar 16 '17 at 19:53

The calculation of $$\int_{-1}^{1} L_k(x)\,\text{d} x$$ for the Lagrange polynomials $$L_k$$ defined on an arbitrary grid $$x_k, k=0,\ldots,n$$ can be performed by the following two steps:

1. Calculate the Clenshaw-Curtis quadrature weights $$w^{\text{cc}}_k$$ on the Chebyshev extrema grid $$y_k$$ for $$k=0,\ldots,n$$: $$y_k = \cos\left(\frac{k\pi}{n}\right)\\ w^{\text{cc}}_k = \frac{c_k}{n }\Bigg(1-\sum_{j=1}^{\lfloor{n/2}\rfloor} \frac{b_j}{4j^2-1} \cos\left(\frac{2\pi \,j \, k}{n}\right)\Bigg)$$ where $$b_k := 1$$ for $$k=n/2$$, otherwise $$b_k:=2$$, and $$c_k:=1$$ for $$k=0$$ or $$k=n$$, otherwise $$c_k:=2$$. See the paper by Waldvogel (2006) for further details.

2. Transform the weights $$w^{\text{cc}}_k$$ to the arbitrary grid $$x_k, k = 0,\ldots, n$$, via the transformation matrix $$M$$ to obtain the sought weights $$w_k$$, $$w_k = \sum_{j} M_{kj} w_j^{\text{cc}}$$ where $$M_{jk} \ = \ L_j(y_k)\,.$$

In principle this is just Clenshaw-Curtis quadrature with function values on the arbitrary grid $$x_k$$, but obtained by basis transformation (for a general refernce on Clenshaw-Curtis, see e.g. the Trefethen paper).

The algorithm seems to be quite stable, particularly when compared to the Vandermonde approach as provided in the answer by @Kirill: although it follows the same ideas -- generate the quadrature weights in a known basis and then transform to the new grid -- this could have been expected as the transformation in terms of the Vandermonde matrix is usually highly ill-conditioned.

Example: Generation of Legendre-Lobatto quadrature weights

We consider the example of Legendere-Lobatto quadrature rule and compare the accuracy to the monomial approach. As a reference, we use the quadrature weights $$w_k^{\text{Leg}}$$ obtained by the Golub-Welsch algorithm for different $$n$$ and calculate the cumulated error $$\epsilon_n = \sum_{k=1}^n \Big(w_k - w_k^{\text{Leg}}\Big)^2$$ Here is the result: One observes that the Clenshaw-Curtis quadrature weights are perfectly stable throughout the considered range of gridpoints and reproduce the Legendre weights up to machine accuracy ($$\sqrt{\epsilon} \sim 10^{15}$$).

Example: Generation of Newton-Cotes quadrature formulas

We consider the generation of Newton-Cotes quadrature formula on equally-spaced grids. Again, one expects an ill conditioning, as, in short, for polynomial interpolation equally-spaced grids are baaad.

In the following picture, I calculated the absolute sum of the weights $$\sum_i |w_i|/N$$. Up to, say, 50 gridpoints, the result of the monomial and Clenshaw-Curtis approach agree. Thereafter, Clenshaw-Curtis becomes better--for what it's worth. A direct interpretation is, that the equally spaced grid ruins right everything for, say, $$n>10$$. At around $$n=50$$, however, the condition of the Vandermonde matrix strikes back and leads to an even worse result.

This example is due to @NicoSchlömer. I didn't know these rules so far, so I took the abscissas from this implementation and applied both the Vandermonde and the transformed Clenshaw-Curtis approach (where, as above, the Vandermonde approach is using the Björk-Pereyra algorithm).

As suggested in the comment I then calculated the error of integrating a constant function by $$\epsilon = \frac{1}{n}\Bigg|2-\sum_{i=1}^n w_i\Bigg|\,,$$ with the following result: From this picture, the transformed Clenshaw-Curtis approach seems way more efficient than the Vandermonde approach (at least in finite-precision arithmetic). Still, Clenshaw-Curtis breaks down starting from index 7, so other methods should be used.

• Thanks for the interesting answer. I played around with it a little bit, but found the round-off to be significant. For example, the sum of the $w_k$ should always be 2. That's true up to n < 5, but for n == 6 already I'm getting 1.9999949955991916, with and error in 7th significant decimal. – Nico Schlömer Dec 18 '19 at 12:53
• I forgot to mention this is with Gauss-Patterson points, github.com/nschloe/quadpy/blob/master/quadpy/line_segment/…. – Nico Schlömer Dec 18 '19 at 13:13
• @NicoSchlömer: ok, that's interesting, because except for the arbitrary-precision arithmetic (which the julia implementation seems to offer), it's hard to reason how the Vandermonde approach should anyhow be more accurate (as the Vandermonde is among the prime examples for an ill-conditioned matrix). I'll make an edit to the question and consider the Gauss-Patterson quadrature. – davidhigh Dec 19 '19 at 6:41
• Thanks David for looking into all of this. It is interesting indeed! I think it would be a worthwhile addition of this post to include a comparison with ordinary Gauss-Legendre quadrature of the appropriate degree. I would guess it performs like your Clenshaw-Curtis approach. Also, putting the code on github or so and linking it from here will be helpful for anyone looking into this in the future. If you link the top-voted answer, I will make this one the "correct" one because of the interesting insights. – Nico Schlömer Dec 19 '19 at 9:05
• @NicoSchlömer: thank you. What do you mean by comparison with Gauss-Legendre? (because the Gauss-Legendre weights are reproduced to machine accuracy). The comparison between Leg. and CC was done by Trefethen, with the result that the accuracy of CC is often comparable. What would be interesting indeed is to study the performance of different custom grids and compare it to Legendre or Clenshaw-Curtis. Moreover, I linked the top-voted answer. – davidhigh Dec 19 '19 at 10:41

You can evaluate this using the Björck-Pereyra algorithm for solving Vandermonde systems, because you are evaluating $b^\top V^{-1}$ with $b=(2,0,\frac23,0,\frac25,0,\ldots)$, and the algorithm is known to be forward-stable (see Error analysis of the Björck-Pereyra algorithms for solving Vandermonde systems by Nick Higham, http://www.maths.manchester.ac.uk/~higham/narep/narep108.pdf).

Note: it seems that this analysis relies on the property that $0\leq x_1<x_2<\cdots <x_n$ (which is equivalent to $V$ being totally positive), as well as the elements of $b$ having alternating signs (which ensures there's no catastrophic cancellation in the subtractions below), in which case the errors are independent of the condition number, and it will still work in the more general case, without $0\leq x_1$, but the error bounds will be different. In any case, it takes $O(n^2)$ time, and avoids the problem with evaluating/integrating $L_i$'s, so it might be worth it even then, but I hadn't realized this point when I started writing this answer. You might be able to just map $x\mapsto \frac12(x+1)$ if that works for your problem.

I wrote a small Julia program to check that this actually works, and it gives $O(\epsilon_{\mathrm{mach}})$ relative errors.

module VandermondeInverse

using SpecialMatrices

function main(n=8)
X = Rational{BigInt}[k//(n-1) for k=0:n-1]
# X = convert(Vector{Rational{BigInt}}, linspace(-1, 1, n))
x = convert(Vector{Float64}, X)

A = convert(Matrix{Rational{BigInt}}, Vandermonde(X))
b = [i%2==0 ? 2//(i+1) : 0 for i=0:n-1]
println("Norm: ", norm(A, Inf))
println("Norm of inverse: ", norm(inv(A), Inf))
println("Condition number: ", cond(convert(Matrix{Float64}, A)))
ans = A'\b

B = convert(Matrix{Float64}, A)
c = convert(Vector{Float64}, b)

println("Linear solve: ", norm((B'\c - ans)./ans, Inf))

d = vec(c')
for k=1:n, l=n:-1:k+1
d[l] -= x[k]*d[l-1]
end

for k=n-1:-1:1, l=k:n
if l > k
d[l] /= x[l]-x[l-k]
end
if l < n
d[l] -= d[l+1]/(x[l+1] - x[l-k+1])
end
end
println("Neville elimination: ", norm((d-ans)./ans, Inf))

nothing
end

end

V = VandermondeInverse


Output:

julia> V.main(14)
Norm: 14.0
Norm of inverse: 1.4285962612120493e10
Condition number: 5.2214922998851654e10
Linear solve: 1.5714609387747318e-8
Neville elimination: 1.3238218572356314e-15


If X isn't positive like in this test, then it seems the relative errors are of the same order as with a regular linear solve.

Why $b^\top V^{-1}$? It's actually a very useful common trick for working with polynomials of all types, but especially the Lagrange interpolating polynomials, converting the problem to a matrix form. The condition that defines $L_i$ is that $L_i(x_j)=\delta_{ij}$. Let $\alpha_{jk}$ be the coefficients of $L_k$, i.e., $$L_k(x) = \sum_{j,k} \alpha_{j,k}x^j = (1,x,x^2,\ldots,x^n)^\top (\alpha_{0k},\ldots,\alpha_{nk}),$$ and $L$ be the whole matrix of coefficients, arranged by columns: $$L = \begin{pmatrix} \alpha_{00}& \cdots & \alpha_{0n}\\\vdots &&\vdots\\ \alpha_{n0}&\cdots&\alpha_{nn} \end{pmatrix}.$$ Because of the definition of $L_k$ above as a vector product, multiplying $L$ on the left by $(1,x,\ldots,x^n)$ yields $$(1,x,x^2,\ldots,x^n)L = (L_0(x),L_1(x),\ldots,L_n(x)).$$ Using the condition $L_k(x_j)=\delta_{jk}$, this means that $$\begin{pmatrix} 1&x_0&x_0^2&\cdots&x_0^n\\ \vdots\\ 1&x_n&x_n^2&\cdots&x_n^n \end{pmatrix} L = I,$$ so $L = V^{-1}$, where $V$ is the Vandermonde matrix of $x_j$'s.

Finally, since $\int_{-1}^{1}x^k\,\mathrm{d}x = \frac{1+(-1)^k}{k+1}$, we have $$\int_{-1}^{1}L_k(x)\,\mathrm{d}x = \sum_j \alpha_{jk}\frac{1+(-1)^k}{k+1} = (2,0,\tfrac23,0,\tfrac25,0,\ldots)^\top (\alpha_{0k},\ldots,\alpha_{nk}).$$ So the $n+1$ numbers you are looking for, for $k=0\ldots n$, are given by $(2,0,\tfrac23,0,\ldots)^\top L$, where $L=V^{-1}$ is the inverse of the Vandermonde matrix.

• Thanks for the elaborate answer! To add even more clarity, could you give more detail on how to bring the problem into Vandermonde form? – Nico Schlömer Mar 18 '17 at 10:26
• @NicoSchlömer Sure, see edit. Thank you for the question, I had no idea this algorithm even existed. – Kirill Mar 18 '17 at 18:06
• I played around with this a little bit. My computer's memory is exhausted by Julia's rationals at $n=16$ where all still is fine. I have one more reference solution for $n=31$, and unfortunately even Björck-Pereyra breaks down then. Pity! It's quite curious that it provides such accurate results in the first place. In any case, thanks for the valuable input. – Nico Schlömer Mar 20 '17 at 19:44
• @NicoSchlömer Hm, I can't reproduce that: with the above code, V.main(32) produces a sensible answer in about a second on my laptop (while using only a little memory). The numbers aren't even that big, the largest numerator has 54 digits, so I suspect something else is going wrong for you. Can you post a gist, because I'm curious to see how it fails? – Kirill Mar 20 '17 at 22:31
• @NicoSchlömer If you mean the way it starts using BigFloat's when printing the relative errors for the output, I'm not quite sure what the rule is. But I made sure it uses Float64 for d: check with @show typeof(d). Let me know if you find more problems with it. – Kirill Mar 20 '17 at 22:50

Calculate the products of the nominators and the denominators first and then divide once. The two products should be of the same order of magnitude, so there should be no significant round-off errors. Also you get the added benefit of increased speed, due to the reduced number of floating point calculations.