# 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

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.