In a problem I working on, I have the following value $$ y = f(a_1,a_2,\ldots,a_n) = \varphi^{a_1} + \varphi^{a_2} + \cdots + \varphi^{a_n} \enspace, $$ for $(a_1,a_2,\ldots,a_n)\in\{0,1,\ldots\}$ and $\varphi = 1 - 1/m$ for $m$ an integer greater than 1. I am supposed, given $y$, to recover $$ x = a_1 + a_2 + \cdots + a_n \enspace.$$ I am able to show that the $x$ is unique when the $a_i$'s are whole numbers. For this question, I consider $x$ to be the "inverse function" of $y$, which makes sense somewhat since they are both symmetric in the $a_i$'s.
I have two solutions but none of them is satisfatory (one is an approximation and the other is computationally demanding and suffers from numerical issues).
Solution 1 (using truncated Taylor expansion): $$ x = (y-n)/\ln(\varphi) + O(m^{-1}) \enspace. $$ The problem with this solution is that the error is high when $m$ is small, and I do not control $m$.
Solution 2 (Tracing the hypersurface):
Basically trying out all the integers lying (close) to the hypersurface defined by the implicit equation of $y$ until we find an exact match (lying exactly on the hypersurface). The problem with this solution is that it can take too much time if the $a_i$'s are huge, even if we exploit the symmetry that the order of $a_i$ doesn't matter. Additionally, $y$ lies between 0 and 2, and tiny changes in $y$ may lead to big changes in $x$ (the function may be numerically unstable). Furthermore, due to the use of floating point representation, we output the point closest to the hypersufrace instead.
One more problem In addition to the above challenges, what I am actually given is a perturbed value of $y$: $y+\text{noise}$. So I am also interested in understanding the behaviour of a given algorithm for extracting $x$ under perturbation.
So my questions are:
Is this function known under a name? Has it (or its "inverse") been studied before? Is there a direct formula (not necessarily "closed-form") to compute the "inverse" than the two ways I proposed? What is the behaviour of the inverse under perturbation?
Here's a conceptual proof that the answer can't be unique, along with a specific example to illustrate. I'll rewrite your (multi)set $\{a_1, a_2, \ldots, a_n\}$ as $\{m_1\cdot d_1, m_2\cdot d_2, \ldots, m_t\cdot d_t\}$; that is, the values $d_i$ (degree) with multiplicity $m_i$ (e.g., the multiset $\{0, 0, 1, 2, 2, 2\}$ with $n=6$ corresponds to $\{2\cdot 0, 1\cdot 1, 3\cdot 2\}$.
Now, consider the polynomial $p(x) = \sum_{i=1}^nx^{a_i}$ $= \sum_{j=1}^tm_jx^{d_j}$.
You say that you know the value $n$; this value is just $p(1)$. You also know the value $p(\varphi)$; and you want to know the value $\sum_i a_i=\sum_j m_jd_j$. This value is just $p'(1)$.
But now, consider $r(x)=p(x)+q(x)$, where we've chosen $q(x)$ such that $q(1)=q(\varphi)=0$. This will have the same values at $1$ and $\varphi$, but it won't (generically) have the same derivative at $1$; this means that your hoped-for sum isn't guaranteed to be unique.
For a specific example, consider $\varphi=\dfrac23$; one candidate polynomial is $q(x)=3(x-1)(x-\varphi) = 3x^2-5x+2$. For simplicity, we can take $p(x)=5x$, so $r(x)=3x^2+2$. In other words, the multisets of $a_i$s are $\{1, 1, 1, 1, 1\}$ and $\{0, 0, 2, 2, 2\}$. These have exactly the same number of elements, and in each case the sum $\sum_i\varphi^{a_i}=\dfrac{10}{3}$. But $\sum_ia_i=5$ for the first multiset and $6$ for the second.
OTOH, if it's known that $m_i\lt m$ for all $i$ and an exact value for $p(\varphi)$ is known, then the full multiset can be recovered. This is because the expression acts as a sort of 'base-$m$ expansion': consider multiplying $p(\varphi)$ by $m^{d_1}$, where (for convenience) $d_1$ is the largest degree. Then the term $m_i\varphi^{d_i}$ becomes $m_im^{d_1-d_i}(m-1)^{d_i}$; this will be divisible by $m$ for all $d_i\lt d_1$, and if $d_i=d_1$ then the expression is congruent to $(-1)^{d_1}m_1\pmod m$; this $p(\varphi)$ in lowest terms is guaranteed to have a denominator of $m^{d_1}$; the expression $m^{d_1}p(\varphi)$ is manifestly integral, and since it's nonzero mod $m$ no smaller power of $m$ will produce an integer. This means that the sequence of $m_i$ and $d_i$ can be extracted by finding the first $d$ with $m^dp(\varphi)$ integral, setting $d_1=d$ and $m_1=(-1)^dm_dp(\varphi)\bmod m$, and subtracting $m_1\varphi^{d_1}$ from $y=p(\varphi)$, and iterating.
Finally, in the regime where $m\gg m_i$ for all $i$, we can take advantage of the special nature of the value we're after to get an excellent approximate solution. Remember that the sought-after value ($x$ in OP's notation) is just $p'(1)$; but for large $m$, we have $\varphi=1-\frac1m\approx 1$, so we can use the classical approximation that $p'(1)=\dfrac{p(1)-p(1-\epsilon)}{\epsilon}$. In the notation of this problem, this says that $x\approx m(n-y)+O(m^{-1})$ (this is a simpler version of OP's Taylor series approximation). Better yet, we can get good bounds on the error: by the Taylor theorem with remainder, we have $x=m(n-y)+m^{-1}p''(\zeta)$ for some $\varphi\leq\zeta\leq 1$. But note that (writing in terms of the $a_i$ for convenience) since all the $a_i$ are positive, $p''(\zeta)\leq p''(1)$ for all such $\zeta$; what's more, the second derivative of $x^{a_i}$ is $\dfrac{a_i(a_i-1)}{2}x^{a_i-2}\leq \dfrac{a_i^2}{2}x^{a_i-1}\leq\frac12a_i^2$. What's more, we have $\sum_ia_i^2\leq (\sum_ia_i)^2$ (we can do much better than this, but this is fine), so we know that $p''(\zeta)\leq \frac12n^2$ (and in fact, it will generally be much less than this value). Thus, the maximum error in the approximation $x\approx m(n-y)$ is $\frac12n^2m^{-1}$; over the regime listed by OP in comments, this is much smaller than 1 and so the nearest integer to $m(n-y)$ will be the correct answer.