Obtain $\{h_1,h_2,\ldots\}$ from $\{f(0),f(1),f(2),\ldots\}$ with $f(s)=\sum_{i=1}^n \exp(-h_i s)h_i$

I need to obtain $ \{h_1,h_2,\ldots\}$ when given $ \{f(0),f(1),f(2),\ldots\}$ with $ f$ defined as follows $ $ f(s)=\sum_{i=1}^n \exp(-h_i s)h_i$ $

We know that $ h_i>0$ . Is there a way to do this numerically in Mathematica?

As a realistic example, take $ h_i=i^{-p}$ with $ p>1$ , so we can produce $ \{h_i\}$ and $ \{f(s)\}_s$ for testing. I’d be very curious to see if the backward method can be fixed in the snippet below to produce correct values. (without relying on knowledge that $ h_i=i^{-p}$ )

p = 2; n = 10; hi = Table[i^-p, {i, 1, n}]; forward[s_] := Total@fixnan@Table[Exp[-h*s] h, {h, hi}]; fi = Table[forward[s], {s, 0, 20}]; (* 20 is arbitrary,can use more *) backward[fi_] := ConstantArray[1, n]; Print["Error: ", N@Norm[hi - backward[fi]]] (* Error:2.82533,but want small *) 

background this would let one get get spectrum of the Hessian just by knowing the sequence of loss values seen during gradient descent (mathoverflow post).