efficient way to compute lagrange polynomials

For two given vectors, not necessairly the same length, $x=(x_1,\dots,x_N)$ and $s=(s_1,\dots,s_M)$ I would like to compute the following term as efficiently as possible in R. For $j=1, \dots, N$

$$H_j:=\sum_{i=1}^Ml_i(x_j)$$ where $$l_i(x_j)= \prod_{l = 1,l \not= j}^M \frac{x_j-s_l}{s_i-s_l}$$

below is some test code which is very inefficient. Howeve, I dont see how I can speed this up using vectorized style or another smart trick.

x <- rnorm(1000)
s <- rnorm(10,2,0.2)

l <- rep(0,length(s))
output <- rep(0,length(x))
for(j in 1:length(x)){
xj <- x[j]
for(i in 1:length(s)){
s.i.dropped <- s[-i]
l[i] <- prod((xj-s.i.dropped)/(s[i]-s.i.dropped))
}
output[j] <- sum(l)
}
• Questions like this are pretty common on Stack Overflow; I suspect it would be on topic there.
– gung
Aug 9 '16 at 16:35

Why do it yourself when experts have already done it for you? See polynom and polynomF for specific implementations in R.

• thanks for you anwser. I've browsed the two packages. One additional question, do you know a package where it is possible to extract just the coefficient functions $l_i(x)$ in the representation of a lagrange polynomial $f(x) = \sum_{i=1}^Ny_il_i(x)$ ?
– math
Aug 10 '16 at 6:48
• I'm not sure if there's a specific call that will give you this, at least with the libraries I suggested. You could exploit the orthogonality and determine what the coefficients must have been on your own (though this might be too costly for your liking). Aug 10 '16 at 17:06

If you want to know more about efficient numerical procedures to compute Lagrange polynomials, I recommend the following reference.

J.-P. Berrut and L. N. Trefethen, “Barycentric Lagrange Interpolation,” SIAM Rev., vol. 46, no. 3, pp. 501–517, 2004.

The main point of the "barycentric" Lagrange formula is to pre-compute barycentric weights (in $O(n^2$) operations) so that the interpolation can then be subsequently carried out in $O(n)$ operations. The "naive" method you have implemented above takes $O(n^2)$ operations for each interpolation and is hence slower, as you have noted. This formula is also numerically stable near the interpolation points, which is a big bonus.

I am not familiar with R, so I don't know if it already implements this routine. However, it may not be hard to implement yourself.

• Many thanks for the reference! will have a look
– math
Aug 16 '16 at 9:00