# 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).