Lagrange Inversion Theorem Proof Note: throughout this question, I'll be using the following notation convention:
$$f^{(n)}(x)=\frac{d^nf}{dx^n}(x)$$
I was browsing through Wikipedia and even MSE's related questions searching for a proof for the Lagrange Inversion Theorem. I'm letting the function $g$ be $f$'s inverse, and I'm letting $f$ be analytic at $a$ while also maintaining the property that $f'(a) \neq 0$.
$g$ is also assumed to be analytic, therefore suggesting that
$$g(f(w))=\displaystyle\sum_{n=0}^{\infty}g^{(n)}(f(a))\frac{(f(w)-f(a))^n}{n!}=w.$$
Expanding past the zeroth term, we get
$$g(f(w))=a+\displaystyle\sum_{n=1}^{\infty}g^{(n)}(f(a))\frac{(f(w)-f(a))^n}{n!}.$$
Now, letting $f(w)=z$, we have
$$g(z)=a+\displaystyle\sum_{n=1}^{\infty}g^{(n)}(f(a))\frac{(z-f(a))^n}{n!}.$$
Wikipedia states that 
$$g(z)=a+\displaystyle\sum_{n=1}^{\infty}\lim_{w \to a}\Bigg[\frac{d^{n-1}}{dw^{n-1}}\Bigg(\frac{w-a}{f(w)-f(a)}\Bigg)^n\Bigg]\frac{(z-f(a))^n}{n!},$$
which suggests the following:
$$g^{(n)}(f(a))=\lim_{w \to a}\Bigg[\frac{d^{n-1}}{dw^{n-1}}\Bigg(\frac{w-a}{f(w)-f(a)}\Bigg)^n\Bigg].$$
It looks like a residue of some weird function at $a$, but the nature of the function itself and the reason why it shows up from a derivative rather than an integral are really stumping me. I've read a few pdf's related to it, but one gave an answer that was too far from the equation itself for me to draw a connection, and the other was in French (my high school French classes have failed me, needless to say).
This is primarily just me wanting to find an inverse polylogarithm, a question that has already been asked here on MSE with no avail, but also for me to understand this theorem more precisely. 
If someone could please point out where this residue-ish thing comes from and why it looks the way it does (i.e., no factorial or $2\pi i$ term), I'd be greatly appreciative. If instead you were able to give me a derivation or a proof of the polynomial approximation given by the commenter Simon here, that would be great too. A few more things that may or may not be useful specifically for the polylogarithms are the following equations:
$$z _{s+1}F_s(1,\dots,1;2,\dots,2;z)=\operatorname{Li}_s(z)$$
$$(z-1) _2F_1(1,1;2;1-z)=\log(z)$$
$$\operatorname{Li}_1(z)=-\log(1-z)$$
$$\operatorname{Li}_0(z)=\frac{z}{1-z}$$
$$\operatorname{Li}_1^{-1}(z)=1-e^{-z}$$
$$\operatorname{Li}_0^{-1}(z)=\frac{z}{1+z}$$
were $_pF_q(a_1,\dots,a_p;b_1,\dots,b_q;z)$ denotes the generalized hypergeometric function and $\log$ denotes the natural logarithm.
 A: I will write what I have found here because it seems interesting enough. It is more a demonstration of how to group the terms in the messy reversion into references to documented (but complicated) sequences.
For the Polylogarithm we have the series representation
$$
\mathrm{Li}_s(z) = \sum_{k=1}^\infty \frac{z^k}{k^s}
$$
if we perform a series reversion on this (term by term) we end up with an expansion for the inverse function
$$
\mathrm{Li}^{-1}_s(z) = \sum_{k=1}^\infty a_k z^k
$$
the first few coefficients are
$$
a_1 = 1 \\
a_2 = -2^{-s} \\
a_3 = 2^{1-2s} - 3^{-s} \\
a_4 = 5 6^{-s} - 8^{-s}(5+2^s) \\
\cdots
$$
there may be a pattern in there somewhere, but the terms seem to grow quite large and complicated rather quickly. For some reason I considered looking at the inverse Mellin transform of these coefficients, multipled by a gamma function, we can denote these functions $e_k(x)$
$$
e_k(x) = \mathcal{M}^{-1}[\Gamma[s]a_k(s)](s)
$$
these begin \begin{equation}
e_1(x) = e^{-x} \\
e_2(x) = -e^{-2x} \\
e_3(x) = 2 e^{-4x} - e^{-3 x} \\
e_4(x) = -5 e^{-8 x} + 5 e^{-6x} - e^{-4 x} \\
\cdots
\end{equation}
in each term $e_k(x)$ there are $P(k-1)$ exponential functions, where $P(k)$ from $k=0$ goes like $1,1,2,3,5,7,\cdots$ and are the partition numbers A000041. The coefficients in this grid of exponentials goes like:
$$
\alpha_1=[ +1]\\
\alpha_2=[ -1]\\
\alpha_3=[ -1, 2]\\
\alpha_4=[ -1, 5, -5]\\
\alpha_5=[ -1, 6,  3, -21,  14]\\
\alpha_6=[ -1, 7,  7, -28, -28, 84, -42]\\
\alpha_7=[ -1, 8,  8,   4, -36, -72, -12, 120, 180, -330, 132]
$$
, and appear to be given by A111785. Interestingly, the powers of the exponentials also appear to have a sequence, the terms go like A074139, which is titled "Number of divisors of A036035(n).". A036035 is titled "Least integer of each prime signature, in graded (reflected or not) colexicographic order of exponents."
We can recreate a coefficient by performing the Mellin transform and dividing through by $\Gamma(s)$
$$
a_k(s) = \frac{1}{\Gamma(s)}\mathcal{M}[e_k(x)](s) = \frac{1}{\Gamma(s)}\int_0^\infty x^{s-1}e_k(x) \; dx
$$
Then we can write
$$
e_k(x) = \sum_{l=1}^{P(k-1)} \alpha_{kl}e^{-\beta_{kl} x}
$$
where $\beta_k$ are the anaolgous rows of A074139
$$
\beta_1=[1]\\
\beta_2=[2]\\
\beta_3=[3,  4]\\
\beta_4=[4,  6,  8]\\
\beta_5=[5,  8,  9, 12, 16]\\
\beta_6=[6, 10, 12, 16, 18, 24, 32]\\
$$
we we know then the Mellin transform of $a e^{-bx}=a b^{-s}\Gamma(s)$, and it's as simple as, grouping the terms which are now understood
$$
a_k(s) = \frac{1}{\Gamma(s)}\int_0^\infty x^{s-1}\sum_{l=1}^{P(k-1)} \alpha_{kl}e^{-\beta_{kl} x} \; dx
$$
$$
a_k(s) = \sum_{l=1}^{P(k-1)}\frac{1}{\Gamma(s)}\int_0^\infty x^{s-1} \alpha_{kl}e^{-\beta_{kl} x} \; dx
$$
$$
a_k(s) = \sum_{l=1}^{P(k-1)}\frac{1}{\Gamma(s)}\alpha_{kl}\beta_{kl}^{-s}\Gamma(s)
$$
$$
a_k(s) = \sum_{l=1}^{P(k-1)}\alpha_{kl}\beta_{kl}^{-s}
$$
$$
\mathrm{Li}^{-1}_s(z) = \sum_{k=1}^\infty a_k(s)z^{k} = \sum_{k=1}^\infty \sum_{l=1}^{P(k-1)}\alpha_{kl}\beta_{kl}^{-s}z^k
$$
the limitation here is understanding where the terms in the sequences referenced come from, and finding whether any tractable forms exist for $\alpha$ and $\beta$. We already know that $\beta_{kl} = \sigma_0(\gamma_{kl})$ for the $\gamma$ in A036035, and divisor counting function $\sigma$.
$$
\mathrm{Li}^{-1}_s(z) = \sum_{k=1}^\infty \sum_{l=1}^{P(k-1)}\frac{\alpha_{kl}z^k}{\sigma^s(\gamma_{kl})}
$$
A: I haven't derived a general case of the theorem, but here's proof for the series expansion centered at $a=0$. Credit: https://users.math.msu.edu/users/magyarp/math880/Lagrange.pdf
$\mathbf{Lagrange\,\, Inversion\,\, theorem}$: Let $f:A\rightarrow B$ be holomorphic in a neighbourhood of $z=0$, and suppose that $f(0)=0$ and $f'(0)\neq 0$ (this is for the Inversion function theorem). Let $C$ be the circle $\partial D(0,\epsilon)$, the circle centered at $0$ with $\epsilon$ radius.
Now let $g:B\rightarrow A$ be the inverse function of $f$, such that $g(f(z))=z$.  Then:
$$g(z)=\sum_{n=1}^{\infty}\frac{c_n}{n!} z^n,\,\,\mathrm{where}\,\, c_n=\lim_{z\rightarrow 0}\frac{d^{n-1}}{dz^{n-1}}\left[\left(\frac{z}{f(z)}\right)^n\right].$$
$\mathbf{Proof}$: Considering $z$ lies inside $f(C)$, Cauchy integral formula gives: $$g(z)=\frac{1}{2\pi i}\oint_{f(C)}\frac{g(\zeta)}{\zeta-z}\, d\zeta=\frac{1}{2\pi i}\oint_C \frac{uf'(u)}{f(u)-z}\, du$$
Now focus on the integrand: $$\frac{uf'(u)}{f(u)-z}=\frac{uf'(u)}{f(u)}\frac{1}{1-\frac{z}{f(u)}}=\frac{uf'(u)}{f(u)}\sum_{n\geq 0}\left(\frac{z}{f(u)}\right)^n, |z|<|f(u)|$$
Notice that the requirement $|z|<|f(u)|$ can be satisfied, since a circle $C$ can be found such that $f(u)\neq 0$.  Since the function $f$ is holomorphic, there must be no non-isolated zeros. Thus we can set $\epsilon$ to be strictly less that the distance of $0$ and all the other zeros of $f(C)$
Putting it back in the integral:
$$g(z)=\frac{1}{2\pi i}\oint_C\frac{uf'(u)}{f(u)}\sum_{n\geq 0}\left(\frac{z}{f(u)}\right)^n\, du=\frac{1}{2\pi i}\sum_{n\geq 0}z^n \oint_C\frac{uf'(u)}{f(u)^{n+1}} du$$
The shape of the power series can already be seen: the coefficient is the contour integral right-hand-side divided by $2\pi i$.  Doing integration by parts (differentiating $u$ and integrating $f'/f^n$):
$$\frac{1}{2\pi i}\sum_{n\geq 0}z^n \oint_C\frac{uf'(u)}{f(u)^{n+1}}\, du=\frac{1}{2\pi i}\sum_{n\geq 0}z^n \oint_C\frac{1}{n f(u)^n}\, du=\sum_{n\geq 0}b_n z^n$$
Where$$b_n=\frac{1}{2\pi i} \oint_C\frac{1}{n f(u)^n}\, du$$ Using the residue theorem, $$b_n=\frac{1}{n}\mathbf{Res} \left(1/f(u)^n,0\right)$$
Since $f(u)$ is has no constant term, it's power series would be like $$f(z)=a_1 z+a_2 z^2+a_3z^3...$$
$$\frac{1}{f(z)^n}=\frac{1}{z^n(a_1+a_2z+a_3z^2...)^n}$$ has a $n$-th order pole.
$$\mathbf{Res} \left(1/f(u)^n,0\right)=\frac{1}{(n-1)!}\lim_{u\rightarrow 0}\frac{d^{n-1}}{dz^{n-1}}\frac{u^n}{f(u)^n}$$
Before putting it back to $b_n$. Note that $b_0$ is zero because by assumption $f(0)=0$, and $g$ is the inverse of $f$ so $g(0)=g(f(0))=0$
$$b_n=\frac{1}{n!}\lim_{u\rightarrow 0}\frac{d^{n-1}}{dz^{n-1}}\frac{u^n}{f(u)^n}$$
Finally putting $b_n$ back into $g(z)$, $$g(z)=\sum_{n\geq 0}b_n z^n=\sum_{n\geq 1}\frac{z^n}{n!}\lim_{u\rightarrow 0}\frac{d^{n-1}}{dz^{n-1}}\frac{u^n}{f(u)^n}$$
