Looking up Lagrange interpolation on Wikipedia, I found something new to me: the barycentric form of Lagrange interpolation. I thought it would be instructional to implement the various forms in q. As usual, my code is available: lagrange_q.
Simple Lagrange interpolation
For simple Lagrange interpolation, the main task is to evaluate the polynomials:
.
The tricky part is the constraint in the product. This could be done using:
q)dij:{x[(c,c-1)#reverse til c:count x]} q)show dij til 5 4 3 2 1 0 4 3 2 1 0 4 3 2 1 0 4 3 2 1 0
I would have preferred
1 2 3 4 0 2 3 4 0 1 3 4 0 1 2 4 0 1 2 3
but “dij” suffices as the product does not care about the exact ordering of the terms. So here is my simple Lagrange interpolation in q:
sili:{[x;y] {[x;y;w;n] sum y*prd flip(n-w)%x-w}[x;y;w:dij x]}
What this function does is to take in 2 lists x and y = F x corresponding to known points of a function F and it produces a polynomial function in n that interpolates F in those known points. For example,
q)F:{1f-2*x*3+4*x*5-6*x} q)xs:2.0 3.5 7.11 13.17 q)ys:F xs q)s:sili[xs;ys]
In the above, s is the polynomial function interpolating the 4 known points of F. Since F is of degree 3, the result is exact (barring small errors in floating point operations).
Modified Lagrange interpolation
The modified Lagrange interpolation rewrites the simple form into:
where
The main efficiency gain over simple Lagrange interpolation is the pre-computation of the weights . The modified form is implemented in q as:
moli:{[x;y] w:1%prd flip x-dij x; {[x;y;w;n] $[n in x;y x?n;prd p,sum y*w%p:n-x]}[x;y;w]}
Due to the possibility of division by 0 when n is an existing point in x, it is necessary to check for “n in x”.
Barycentric Lagrange interpolation
The barycentric form of Lagrange interpolation goes a step further and removes using the observation that the constant function 1 satisfies
.
The barycentric form of Lagrange interpolation is thus:
The symmetry in this formulation makes it easy to implement in q:
bali:{[x;y] w:1%prd flip x-dij x; {[x;y;w;n] $[n in x;y x?n;sum[w*y]%sum w%:n-x]}[x;y;w]}
Performance
I compared the performance of all 3 forms of Lagrange interpolation:
q)s:sili[xs;ys]; m:moli[xs;ys]; b:bali[xs;ys] q)domain:"f"$til 100000 q)\t s each domain 733 q)\t m each domain 250 q)\t b each domain 221
Much of the performance gain is due to the pre-computation of the weights.
Numerical stability
It is a well-known fact (at least to numerical analysts, of which I am not one) that the choice of interpolating points is important to numerical stability. To illustrate this fact, I use the function . Applying Lagrange interpolation with 97 equally spaced points in the interval [-0.5,0.5], I obtained:
n F Es Em Eb -------------------------------------------------------------- -1 -0.5 7.634519e+63 7.634519e+63 2.873049e+12 -0.8947368 -0.3531856 5.032181e+58 5.032181e+58 3.921853e+12 -0.7894737 -0.2285319 4.53236e+52 4.53236e+52 3.274189e+12 -0.6842105 -0.1260388 2.04947e+45 2.04947e+45 3.413191e+12 -0.5789474 -0.04570637 4.09622e+35 4.09622e+35 3.463492e+12 -0.4736842 0.01246537 8.099713e+18 8.099713e+18 3.103591e+12 -0.3684211 0.04847645 1.360591e+08 1.360591e+08 1.360531e+08 -0.2631579 0.06232687 49.44574 49.44574 49.44574 -0.1578947 0.05401662 0.006490092 0.006490092 0.006490092 -0.05263158 0.02354571 9.073315e-05 9.073315e-05 9.073315e-05 0.05263158 0.07617729 9.073315e-05 9.073315e-05 9.073315e-05 0.1578947 0.2119114 0.006490092 0.006490092 0.006490092 0.2631579 0.3254848 49.44574 49.44574 49.44574 0.3684211 0.4168975 1.360591e+08 1.360591e+08 1.360531e+08 0.4736842 0.4861496 8.099713e+18 8.099713e+18 2.61004e+12 0.5789474 0.533241 4.09622e+35 4.09622e+35 2.911981e+12 0.6842105 0.5581717 2.04947e+45 2.04947e+45 3.075629e+12 0.7894737 0.5609418 4.53236e+52 4.53236e+52 3.078241e+12 0.8947368 0.5415512 5.032181e+58 5.032181e+58 3.082528e+12 1 0.5 7.634519e+63 7.634519e+63 3.694835e+12
The Es, Em and Eb columns show the errors using simi, moli and bali respectively. Notice the growth of the errors near the boundaries of [-0.5, 0.5]. The interpolating polynomial is pretty much useless outside the interval.
Instead of equally spaced points, I also tried Chebyshev points (of the 2nd kind, to be exact). In this case, I obtained:
n F Es Em Eb -------------------------------------------------------------- -1 -0.5 2.97847e+50 2.97847e+50 1.58466e+12 -0.8947368 -0.3531856 9.746906e+44 9.746906e+44 1.444789e+12 -0.7894737 -0.2285319 2.871685e+38 2.871685e+38 1.637349e+12 -0.6842105 -0.1260388 1.719823e+30 1.719823e+30 2.111364e+12 -0.5789474 -0.04570637 2.90491e+18 2.90491e+18 2.650901e+12 -0.4736842 0.01246537 3.764609e-06 3.764609e-06 3.764609e-06 -0.3684211 0.04847645 6.511294e-05 6.511294e-05 6.511294e-05 -0.2631579 0.06232687 2.702486e-05 2.702486e-05 2.702486e-05 -0.1578947 0.05401662 0.000137938 0.000137938 0.000137938 -0.05263158 0.02354571 0.0004967719 0.0004967719 0.0004967719 0.05263158 0.07617729 0.0004967719 0.0004967719 0.0004967719 0.1578947 0.2119114 0.000137938 0.000137938 0.000137938 0.2631579 0.3254848 2.702486e-05 2.702486e-05 2.702486e-05 0.3684211 0.4168975 6.511294e-05 6.511294e-05 6.511294e-05 0.4736842 0.4861496 3.764609e-06 3.764609e-06 3.764609e-06 0.5789474 0.533241 2.90491e+18 2.90491e+18 3.704316e+11 0.6842105 0.5581717 1.719823e+30 1.719823e+30 5.479875e+11 0.7894737 0.5609418 2.871685e+38 2.871685e+38 6.762965e+11 0.8947368 0.5415512 9.746906e+44 9.746906e+44 7.625274e+11 1 0.5 2.97847e+50 2.97847e+50 8.291047e+11
Using Chebyshev points, the interpolating polynomial is stable within the interval [-0.5,0.5].
References
The reader who is interested in the theory of barycentric Lagrange interpolation may consult these excellent papers:
- http://www.comlab.ox.ac.uk/nick.trefethen/barycentric.pdf
- http://www.maths.manchester.ac.uk/~nareports/narep440.pdf