interpolation | thesweeheng's weblog

Lagrange Interpolation in q

31 July, 2008 — thesweeheng

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:

$\displaystyle L(n) = \sum_{j=0}^{k} l_j(n) y_j\text{ where }l_j(n) = \prod_{i=0,i\not=j}^{k} \frac{n-x_i}{x_j-x_i}$ .

The tricky part is the $i\not=j$ 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

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:

$\displaystyle L(n) = l(n) \sum_{j=0}^{k} \frac{w_j}{x-x_j} y_j$

where

$\displaystyle l(n) =\prod_{i=0}^{k} (n-x_i) \text{ and } w_j = \prod_{i=0,i\not=j}^{k} \frac{1}{x_j - x_i}$

The main efficiency gain over simple Lagrange interpolation is the pre-computation of the weights $w_j$ . 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 $l(n)$ using the observation that the constant function 1 satisfies

$\displaystyle 1 = l(n) \sum_{j=0}^{k} \frac{w_j}{x-x_j}$ .

The barycentric form of Lagrange interpolation is thus:

$\displaystyle L(n) = \frac{\sum_{j=0}^{k} \frac{w_j}{n-x_j} y_j}{\sum_{j=0}^{k} \frac{w_j}{n-x_j}}$

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 $F(x) = |x| + x/2 - x^2$ . 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:

thesweeheng’s weblog

Recent Posts