This website shows an implementation in R that fits a Padé approximant:
$$ R(x)= \frac{\sum_{j=0}^m a_j x^j}{1+\sum_{k=1}^n b_k x^k}=\frac{a_0+a_1x+a_2x^2+\cdots+a_mx^m}{1+b_1 x+b_2x^2+\cdots+b_nx^n} $$ to some data.
They write the function like this: $$ y = \frac{a_0 + a_1 x + a_2 x^2}{1 + b_1 x + b_2 x^2} $$
A bit of rearranging:
$$ y (1 + b_1 x + b_2 x^2) = a_0 + a_1 x + a_2 x^2 $$ $$ y = a_0 + a_1 x + a_2 x^2 - b_1 x y - b_2 x^2 y $$
Now, they retrieve the coefficients using R's lm() function (linear model).
Now we have a form that lm can work with. We just need to specify a set of inputs that are powers of x (as in a traditional polynomial fit), and a set of inputs that are y times powers of x. This may seem like a strange thing to do, because we are making a model where we would need to know the value of y in order to predict y. But the trick here is that we will not try to use the fitted model to predict anything; we will just take the coefficients out and rearrange them in a function. The fit_pade function below takes a dataframe with x and y values, fits an lm model, and returns a function of x that uses the coefficents from the model to predict y:
fit_pade <- function(point_data){
fit <- lm(y ~ x + I(x^2) + I(y*x) + I(y*x^2), point_data)
lm_coef <- as.list(coef(fit))
names(lm_coef) <- c("a0", paste0(rep(c('a','b'), each=2), 1:2))
with(lm_coef, function(x)(a0 + a1*x + a2*x^2)/(1 - b1*x - b2*x^2))
}
Some examples are available in the link.
Is there a name for this method of fitting? Is there a paper, article, book, that describes this, either as is, or in more general terms? Preferably something citable?
Bonus question - can this be safely expanded to higher order polynomials?
