Inverse of the polylogarithm The polylogarithm can be defined using the power series
$$
\operatorname{Li}_s(z) = \sum_{k=1}^\infty {z^k \over k^s}.
$$
Contiguous polylogs have the ladder operators
$$
\operatorname{Li}_{s+1}(z) = \int_0^z \frac {\operatorname{Li}_s(t)}{t}\,\mathrm{d}t\,,
\qquad
\operatorname{Li}_{s-1}(z) = z \,{\partial \operatorname{Li}_s(z) \over \partial z}\ ,
$$
and the sequence can be started with either
$$
\operatorname{Li}_{1}(z)  = -\ln(1-z)\,,\qquad
\operatorname{Li}_{0}(z)  = {z \over 1-z} \ .
$$
Both $\operatorname{Li}_0$ and $\operatorname{Li}_1$ have inverse functions (up to a choice of branchcut)
$$
\operatorname{Li}_0^{-1}(z)=\frac{z}{z+1}\,,\quad
\operatorname{Li}_1^{-1}(z)=1-e^{-z}\,,
$$
$$
\operatorname{Li}_0\left(\frac{z}{z+1}\right)
=z=
\operatorname{Li}_1\left(1-e^{-z}\right) + 2 n \pi i\,,\quad n\in\mathbb{Z}
$$
Is there a nice/useful inverse function for the dilog ($\operatorname{Li}_2(z)$) and higher polylogs?
 A: In astrophysics, specifically in partially degenerated matter, are used what is called Fermi-Dirac Integrals, which are written in terms of polylogaritms, and the z-value is a degeneracy parameter. In some papers I found that in fact they need the inverse of the Fermi-Dirac Integrals, that is, the inverse of the Polylogarithm. 
A: I wrote this in a related question and thought I might as well put it here even though the question is quite old. It gives something like a formal representation of the inverse as a series. 
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(\gamma_{kl})^s}
$$
A: I just ran into a need to solve this very problem. I was able to derive a solution using the inversion of a series.
Starting with
$Li_{s}(z) = \sum_{k=1}^{\infty} k^{-s} z^{k}$
we want the inverse function to be of the form
$Li^{-1}_{s} = \sum_{k=1}^{\infty} b_{k} z^{k}$.
Using the inversion of a series technique (see Morse and Feshbach pg 411) I find
$b_k = \frac{1}{k!}\left[\frac{d^{k-1}}{dx^{k-1}}\left(1 + \sum_{n=1}^{\infty} (n+1)^{-s} x^{n}\right)^{-k}\right]_{x=0}$
These coefficients can be fairly easily calculated by hand or using your favorite software suite. As for a deeper pattern, the one interesting thing I found is that the sum in the coefficient is the Lerch transcendent function
$\sum_{n=1}^{\infty} (n+1)^{-s} x^{n} = x \Phi(x,s,2)$
where
$\Phi(x,s,a) = \sum_{n=0}^{\infty}(n+a)^{-s} x^{n}.$
So the coefficients of the inverse polylog are
$b_k = \frac{1}{k!}\left[\frac{d^{k-1}}{dx^{k-1}}\left(1 + x \Phi(x,s,2)\right)^{-k}\right]_{x=0}$
At this point I just calculated the coefficients using Mathematica and got the following for the first few coefficients
$b_1 = 1$
$b_2 = -2^{-s}$
$b_3 = 2^{1-2s} - 3^{-s}$
$b_4 = -24^{-s} (5\times(3^{s}-4^{s}) + 6^{s})$
They quickly get more complex from there.
A: The answer is yes.  Polylogarithms are positive for all parameter values s and all z > 0.  Thus the integrals used in the ladder operators should all have well defined inverses when the argument is positive, as the integrands are all monotonically increasing functions of z.  It strikes me as similar to how the Elliptic Functions are also defined as inverses of integrals of certain functions.  Analytic continuation should then extend this to the whole complex plane (assuming a branch cut somewhere along the negative real axis), though as @Jeffrey-Hersh points out using symbolic series inversion techniques may be complicated.  I wonder if there are other functional equations we could exploit?
How to numerically compute this inverse generically even for positive z seems like another question entirely.  Its likely its own class of special function.  That said, numerically it shouldn't be that hard to implement a low order inverse like $Li_2^{-1}(x)$.  If you have numerical versions of polylog that can handle fractional orders, the method of numerical integral inversion should work for those too, assuming fractional polylogs are also monotonically increasing on the relevant domains.
I fear you'll probably have to write your own functions in whatever symbolic computing environment you're using, as I'm unaware of any libraries specialized to this problem.  If you find an implementation, please share!
A: Answering this for complex $z$ on the open unit disk, for $s>0,$ the polylogarithm is a convex function which is a special subclass of univalent (one to one) functions. It has an inverse but such an inverse (as shown above) is not analytically tractable.
Source: https://academic.oup.com/jlms/article-abstract/s2-27/3/435/828376
