1

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)
}
math
  • 101
  • 1
  • 13
  • Questions like this are pretty common on [SO]; I suspect it would be on topic there. –  Aug 09 '16 at 16:35

2 Answers2

2

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 06: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
2

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.

Nick C.
  • 188
  • 5