# What algorithm is used in Matlab's pade function?

I try to find out for quite some time now, how Matlab implements the calculation of Padé Approximants using its symbolic pade function. (the code of is buried in a compiled mex-file)

I compared its output with the "direct calculation" via coefficient comparison and solving the resulting linear systems of equations. I used floating point, variable precision arithmetic and pure symbolic variables and even in the last case Matlab's implementation achieves far better results for orders $> 7$.

So there must be a more sophisticated algorithm involved, but the documentation remains silent. The only tiny bit of information I found was the "Euclidian Algorithm" - but in my understanding that does not directly lead to the Padé-Approximant, just to the greatest common divisor of two numbers, am I right?

I contacted the Matlab Support and received the answer, that if no source is given the algorithm would be their own development. I actually can't believe that they developed an algorithm from scratch, if there are books full of alternatives (e.g. Baker, see below), or do they? At least their algorithm should be based on something?

Baker, George A.jun.; Graves-Morris, Peter, Padé approximants., Encyclopedia of Mathematics and Its Applications. 59. Cambridge: Cambridge Univ. Press. xiv, 746 p. (1996). ZBL0923.41001.

Can you give me some hints on how they could have implemented the algorithm? What approach for calculating Padé approximants could be based on the mentioned "Euclidian Algorithm"? How could I proceed from there?

I actually could just use their implementation, but I feel for any publication one should know how it is done.

• You can check out the documentation of [expm]. Matlab's code for matrix exponential is based on padé-approximation. There are also expmdemo-files that should show how padé is used. Maybe there are some hints in there? : ch.mathworks.com/help/matlab/ref/expm.html – mdot Dec 31 '17 at 10:40
• Did you check out en.wikipedia.org/wiki/Pad%C3%A9_approximant? There I explained once how to use the extended Euclidean algorithm to get a diagonal of the Padé table. I think I wrote it on the talk page and someone else transferred it with proper references to the article page. – LutzL Dec 31 '17 at 10:47
• @max I already read about the implementation of expm, I will have a deeper look in the algortihm. But I guess the case of $e^A$ is much simpler than the generic case, probably the reason why Matlab does not use a real Padé approximation within the code (it does not call its own pade function), but rather hard coded (!!!) Padé coefficients. – thewaywewalk Dec 31 '17 at 11:13

Citing myself on https://en.wikipedia.org/wiki/Talk:Pad%C3%A9_approximant, one of the classical ways to compute Padé approximants uses the extended Euclidean algorithm. To compute the $[m/n]$ approximant, set \begin{align} &&r_0&=x^{m+n+1},&v_0&=0,\\ &&r_1&=f(x), &v_1&=1,\\ q_1&=r_0\text{ div }r_1,& r_2&=r_0-q_1r_1,& v_2&=v_0-q_1v_1\\ &\ \ \vdots\\ q_k&=r_{k-1}\text{ div }r_k,&r_{k+1}&=r_{k-1}-q_kr_k,&v_{k+1}&=v_{k-1}-q_kv_k, \\ &\ \ \vdots \end{align} where $f(x)$ is the Taylor polynomial of degree $m+n$. Then in every step $$v_kf(x)=r_k\mod x^{m+n+1}$$ which means that $\frac{r_k}{v_k}$ is the $[\deg r_k/\deg v_k]$ Padé approximant. For most $f$ the quotients $q_k$ will all have degree $1$ so that $\deg v_k=k-1$ and $\deg r_k=m+n+1-k$ so that in step $k=n+1$ one obtains the $[m/n]$ approximant.