So next in my list of (my) useless ideas here's the idea of an "inverse Taylor series"...enjoy!
Let $f$ be an analytic function, so that \begin{equation} f(x)=\sum_{n=0}^\infty a_n x^n \end{equation} The idea is that since \begin{equation} f(kx)=\sum_{n=0}^\infty a_n k^nx^n=:f_k(x) \end{equation} one could try to rewrite a power of $x$, say $x^m$, in terms of the $f_k(x)$'s: \begin{equation} \sum_{k\in\mathbb{Z}} c_k f_k(x) = x^m \end{equation} Noting that \begin{equation} \sum_{k\in\mathbb{Z}} c_k f_k(x) = \sum_{k\in\mathbb{Z}} c_k \sum_{n=0}^\infty a_nk^nx^n = \sum_{n=0}^\infty \left(\sum_{k\in\mathbb{Z}} c_ka_nk^n\right)x^n \end{equation} we would have the condition \begin{equation} \sum_{k\in\mathbb{Z}} c_k f_k(x) = x^m \qquad\Longrightarrow\qquad \sum_{k\in\mathbb{Z}} c_ka_nk^n = \delta_{mn} \end{equation} which is an infinite-dimensional linear system in the unknowns $c_k$. It reminds me in some weird way of a Fourier expansion.
As an example, I tried to "reverse-expand" $f(x)=x$ with a series of exponentials (truncating their Taylor expansion at the fourth power of $x$ because I don't remember how to handle infinite-dimensional linear systems).
$$e^{-2x}=1-2x+2x^2-\frac{4}{3}x^3-\frac{2}{3}x^4$$ $$e^{-x}=1-x+\frac{1}{2}x^2-\frac{1}{6}x^3+\frac{1}{24}x^4$$ $$e^{0}=1\qquad\text{(duh)}$$ $$e^x=1+x+\frac{1}{2}x^2+\frac{1}{6}x^3+\frac{1}{24}x^4$$ $$e^{2x}=1+2x+2x^2+\frac{4}{3}x^3+\frac{2}{3}x^4$$
so we have \begin{equation} x=c_{-2}e^{-2x}+c_{-1}e^{-x}+c_0e^0+c_1e^x+c_2e^{2x}= (c_{-2}+c_{-1}+c_0+c_1+c_2)+(-2c_{-2}-c_{-1}+c_1+2c_2)x+\left(2c_{-2}+\frac{c_{-1}}{2}+\frac{c_1}{2}+2c_2\right)x^2+\left(-\frac{4}{3}c_{-2}-\frac{c_{-1}}{6}+\frac{c_1}{6}+\frac{4}{3}c_2\right)x^3+\left(-\frac{2}{3}c_{-2}+\frac{c_{-1}}{24}+\frac{c_1}{24}+\frac{2}{3}c_2\right)x^4 \end{equation}
This yields a $5x5$ linear system in the unknowns $c_{-2},c_{-1}, c_0, c_1, c_2$ which has a unique solution:
\begin{equation} c_{-2}=-\frac{1}{4},\qquad c_{-1}=\frac{2}{3}, \qquad c_0=-2, \qquad c_1 = 2, \qquad c_2 = -\frac{5}{12} \end{equation} The resulting function is \begin{equation} f(x)=-\frac{1}{4}e^{-2x}+\frac{2}{3}e^{-x}-2+2e^x-\frac{5}{12}e^{2x} \end{equation} graphed here along with the desired function $x$:
which "feels plausible" since we've worked with the truncated series...
My question(s):
(1) How do I approach infinite-dimensional linear systems of this type and are there any convergence problems or other dangers in naively treating these things?
(2) Does this silly idea already exist in literature?
EDIT
As a complement to the question here's the graph of the approximation to $f(x)=x$ given by Claude Leibovici's method, which yields
\begin{equation}
g(x) = \frac{1}{12}e^{-2x}-\frac{2}{3}e^{-x}+\frac{2}{3}e^{x}-\frac{1}{12}e^{2x} = \frac{1}{3}\sinh(x)(4-\cosh(x))
\end{equation}
Suppose that you write $$x=\sum_{k=-n}^{k=+n} c_k\, e^{k x}$$ instead of Taylor expansions, consider the norm $$\Phi_n=\int_{-a}^{+a}\Bigg[x-\sum_{k=-n}^{k=+n} c_k\, e^{k x}\Bigg]^2\,dx$$ and minimize it with respect to the $c_k$'s. If we use a symmetric interval, we shall have $c_{k}=-c_{-k}$ and $c_0=0$. So, for the example you gave, consider $a=1$ and $n=2$ and the minimization will lead to $$c_{-2}=-\frac{3 e^2 \left(7-84 e^2+54 e^4-28 e^6+3 e^8\right)}{1+84 e^2-57 e^4+448 e^6-201 e^8+12 e^{10}+e^{12}}$$ $$c_{-1}=\frac{6 e \left(3+9 e^2+38 e^4+6 e^6-9 e^8+e^{10}\right)}{1+84 e^2-57 e^4+448 e^6-201 e^8+12 e^{10}+e^{12}}$$ $$\Phi_2=\frac{3 \left(15+324 e^2+134 e^4+804 e^6-405 e^8+40 e^{10}\right)}{1+84 e^2-57 e^4+448 e^6-201 e^8+12 e^{10}+e^{12}}-\frac{23}{6}=2.75\times 10^{-6}$$ As usual, the maximum error is at the bounds; in this case, it is $3.83\times 10^{-3}$.
Using $n=3$, the norm decreases to $7.37\times 10^{-9}$ (a factor of $373$)
Edit
Using as you did the series expansion $$e^{k x}=1+k x+\frac{k^2 x^2}{2}+\frac{k^3 x^3}{6}+\frac{k^4 x^4}{24}+O\left(x^5\right)$$ and using the same norm to minimize the error, I find $$x=\frac 1{12}e^{-2x}-\frac 23 e^{-x}+\frac 23 e^{x}-\frac 1{12}e^{2x}=\frac{1}{3} \sinh (x) (4-\cosh (x))$$
Here is a plot with the two approximations (black: minimizer of the norm with exponential functions; red: minimizer of the norm using polynomial expansion of $e^{kx}$)
Update
Taking into account the above, it means that we make the approximation $$x_{(n)}=\sum_{k=1}^n c_k \, \sinh(k x)$$ and the parameters are obtained by the minimization of the norm $$\Phi_n=\int_{-a}^{+a}\Bigg[x-\sum_{k=1}^n c_k \, \sinh(k x)\Bigg]^2\,dx$$ where the hyperbolic sines are either used per se or approximated by their Taylor series; the last solution is obviously simpler in terms of apparent "simplicity" of the coefficients.
However, for a given $n$, the results will depend on the number of terms used in the Taylor series. Expanding the hyperbolic sines to $O\left(x^{2p+1}\right)$, the coefficients are (still for $a=1$) $$\left( \begin{array}{ccc} p & c_1 & c_2 \\ 2 & 1.333333333 & -0.1666666667 \\ 3 & 1.266696757 & -0.1365320137 \\ 4 & 1.260198645 & -0.1337118902 \\ 5 & 1.259839496 & -0.1335590439 \\ \cdots &\cdots &\cdots \\ \infty &1.259826765 &-0.1335536879 \end{array} \right)$$