# How to invert a lagrangian polynomial

I'm reading the following paper (Grezlak and Oosterlee) and I have a specific question to a sentence on page 5. I quote:

"Since the mapping $y=g(x)$ is bijective and $g(x)$ is strictly increasing, so is $g^{-1}(y)$. This implies that the arguments $x$ can be obtained by the inverse interpolation of $g(x)$ against $y$, which can be done at essentially no cost."

For given points $x_i,y_i$ they define $g$ to be the lagrange polynomial:

$$g(x) = \sum_{i=1}^Ny_il_i(x)$$ where $l_i(x)=\prod_{j=1,j\not=i}\frac{x-x_j}{x_i-x_j}$. They cite a paper (E. Michalup. On inverse linear interpolation. Scandinavian Actuarial Journal 1950) about this inverse interpolation which is unfortunatelly not available for free. I'm working with R and was browsing the web about inverse linear interpolation and found the following function. However, it seems to be quite costly. Even worse it just takes a single numeric value as input (yval), which means I need to call the function as many time as inverse I want to calculate. Is there another way to do it or a better function in R?

• In english wikipedia you can read an improve with barycentric. Other alternative is Hermite polynomial. You can use divided differences. Look for packages in R to do that. Aug 12 '16 at 17:38

Without knowing anything about inverse interpolation, this problem should be a straight-forward numerical task given the information about $g(x)$.
Firstly, $g(x)$ is said to be monotonously increasing. Therefore, given $y$, find two nodes s.t. $g(x_i)\leq y\leq g(x_{i+1})$ [this is simple since all is ordered].
From there, use any root-finder which should converge rapidly. This method has the drawback that you have to solve a system for every point $y$. What you can do to speed it up is to add any evaluated point of $g(x)$ in the root-finding process to the list which will decrease the number of future iterations assuming you start from a low-degree polynomial, i.e., small $N$.