1
$\begingroup$

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)
}
$\endgroup$
1
  • $\begingroup$ Questions like this are pretty common on Stack Overflow; I suspect it would be on topic there. $\endgroup$
    – gung
    Aug 9 '16 at 16:35
2
$\begingroup$

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

$\endgroup$
2
  • $\begingroup$ 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)$ ? $\endgroup$
    – math
    Aug 10 '16 at 6:48
  • $\begingroup$ 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). $\endgroup$ Aug 10 '16 at 17:06
1
$\begingroup$

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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.