Integrating Lagrange polynomials with many nodes, round-off

Question

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?

Answer 1

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.

---

Example: Guass-Patterson quadrature

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.

Answer 2

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
  println("True answer: ", ans)

  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
True answer: Rational{Int64}[3202439130233//2916000,-688553801328731//52390800,19139253128382829//261954000,-196146528919726853//785862000,6800579086408939//11642400,-43149880138884259//43659000,32567483200938127//26195400,-7339312362348889//6237000,48767438804485271//58212000,-69618881108680969//157172400,44275410625421677//261954000,-2308743351566483//52390800,11057243346333379//1571724000,-209920276397//404250]
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.

Answer 3

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.