interestingly, fln is the same (up to a constant) as applying RT again to time and squaring it: time2 = Pi/4*(Table[RT[time, i, ns], {i, Length@time}]^2); Norm[time2 - fln] // Chop (*outputs 0*)

Just a comment on my previous comment, I actually was testing this with smaller n (n=100) and realized they aren't perfectly equal when n=400. The relative difference though Norm[time2 - fln]/Norm[fln] is on the order of $10^{-7}$ however so they are almost the same.

@ydd The idea of applying RT to time seems interesting. As if an approximation to the inverse is almost the same formula. What I really want is to find a way to compute fln (fl3n) and compare it to fl (fl3), which contains the "correct" values of $x_j$, for each scale. I should have mentioned, but here $x$ is analogous to the fls and $y$ is time. Nonetheless, your observation could potentially hint at a better fitting procedure than the pure inverse power law I use to compute fln and fl3n. Following this comment, I will replace time and fl with y and x, respectively

So this is quite interesting indeed. If I set x3 = Pi/4*(Table[RT[5*y, i, ns], {i, Length@y}]^2) // Quiet;, I get the following Mean[x2 - x3] (Out[]= -3.88663*10^-22). So x3 "is" x2. I think this is surprising (and promising), but I am not yet sure how it helps in determining the exact x2n corresponding to 5*y, as that is my goal. Any ideas?

I didn't want to put this as answer because it's not a general solution: for large valued x, $y_j \approx \frac{1 - e^{-x_j}}{x_j}$ Solving for $x_j$ yields $x_j \approx \frac{1 + y_j \mathrm{ProductLog}[-\frac{e^\frac{-1}{y_j}}{y_j}]}{y_j}$ you can test this by multiplying x by a large value: inv[yj_] = First@SolveValues[(1 - Exp[-xj])/xj == yj, xj]; xBig = 10^6 x; yBig = Table[RT[xBig, i, ns], {i, Length@xBig}] // Quiet; xnBig = inv[yBig]; Norm[xnBig - xBig]/Norm[xBig] (about 10^-7) Ifyou plot xnBig, you'll see it captures all of the behavior of xBig

The approximation $𝑦_j \approx \frac{1-e^{-x_j}}{x_j}$ only works for large $x_j$ though

@ydd That is quite interesting. I am particularly interested in the dynamics for small $x_j$, do you think any approximation could hold for this case?

I wrote this just in case you need it (but I'm sure you're familar with gradient descent) but this is a gradient descent algorithm for your problem, using $\frac{\pi}{4}y^{-2}$ as an initial condition. You can optimize the step size on each descent step if you like. It's not an analytical solution, so I again did not want to post this as an answer wolframcloud.com/obj/dtrimas/Published/gradDescent.nb

@ydd This is wonderful. I will try it on the 5*y dataset. We could potentially think of ways to optimize it, given that my actual data will have around 1000 points. I agree this is not an analytical approach, but a relatively fast gradient descent implementation would be great! It's a relatively annoying expression because we are very close to an actual analytical answer, but perhaps we could achieve some combination of numerical and theoretical approaches.

For example, ns does not necessarily have to be Floor[n/2]. In fact, it might be the case that it can be much lower, given that 'far away' points do not contribute as much to the value of $y_j$. I am not too sure, however, how this would impact the code and its efficiency.

Ok, I made a new version with 'ns=100` instead of 200 and it is basically the same, but 4 times faster. I also realized I was calculating the sum function fp multiple times for the same input value with the way I defined grad and vec as separate functions, so this is quite a bit faster to run: wolframcloud.com/env/dtrimas/Published/fasterGradDescent.nb

@ydd This looks great. Just a question. Is there any way of avoiding negative values of $x_j$? Perhaps the solution is not unique, but $x_j<0$ would not make much sense in my particular modeling context. Any ideas? Perhaps setting Abs@pl[l_List, i_Integer] where needed?

See here, for example.

After a little bit of reading through the documentation, I was able to get FindRoot to work with the constraint $x_j\geq0$ and I added it to the bottom of this notebook: wolframcloud.com/env/dtrimas/Published/fasterGradDescent.nb It's kind of messy the way I did it and probably can be done in a better way. I tested it with n=50 and ns=25 because it only takes 2 s to run 10 iterations with that. I am running it now with n=400 and with ns=100 to see how long it takes.

Thanks! I will take a look. The BFGS algorithm might also be interesting

It looks my method using FindRoot with constraints actually does not produce a good solution, so I may have to look at that.

@ydd Any updates? BFGS seems somewhat trickier to implement but, in the end, all I really need is the best possible approximation given $x_j>0, \forall j$.

No need to implement BFGS. FindRoot works fine and returns a good solution if you specify a small lower bound that isn't exactly zero (like 10^-10) wolframcloud.com/env/dtrimas/Published/fasterGradDescent.nb Also if you use FindMinimum instead and select "QuasiNewton" method it actually uses BFGS anyways (convergence is quite a bit slower however with QuasiNewton compared to Newton). I didn't use this though because FindRoot worked fine once I changed the lower bound from 0 to 10^-10

@vdd Thanks! Just a quick question, what is val? Also, is there a way to implement the bound imposed in FindRoot into your original method?

'val' is never assigned a value, it is just used as an indexing variable in xIndexed which is the variable list sent into FindRoot (xIndexed looks like {val[1],val[2],val[3],...} and then when FindRoot runs it returns a list of replacement rules {val[1]->number1,val[2]->number2,...} where {number1,number2,...} is the solution. As for adding bounds to my original gradient descent, I am not sure if I can implement this, but one thing you could do is at each step, if the step takes below 0, just clip the value to 0. I don't know if this will effect the quality of the solution however.

Are you able to try the non-negative fitting with y equal to this data? Is it relatively fast?

It took ~40 mins. wolframcloud.com/env/1b2e1cb6-8b95-43ef-b0eb-0a16f68cf504 There are artifacts near the endpoints, but otherwise it is a good fit. My guess is that these artifacts come from the fact that I reduced the size of ns to get this to run in an appreciable amount of time. I can run it with the full sized ns=Floor[n/2] tonight to see (will take ~3hrs). I will let you know. I pasted the resulting reconstructed x, along with the y generated from the reconstructed x at the bottom so you can look at them without having to re-run the code

This looks great! I have one final request if you have the time. The actual data I want to fit is 5*y, where y is linked in my previous comment. Furthermore, there is one slight change to the main function, which might play a crucial role. Change the definition of RT to include the factor $1/v$ on the exponents of the numerator, as seen in this expression, where $v=1.4$. If you could try fitting that, I would really appreciate it, because this "rescaling" might yield an upper bound on y for positive x, set by the list size, and that is an issue.

If that does not converge, try with a 2000-point dataset, given here. I think the issue with this scale, is that RT, as x_j tends to to 0, has a bound that is lower than some dataset values. For example, if you apply fp to ConstantArray[10^(-10), 1000], the result is always lower than some value, which falls lower than values in y, for example.

By increasing the number of points, this bound increases as well. In particular, for 2000 points, the bound is greater than 500, so there is hope that the provided data (which is always less than 500, for example) can be fitted with non-negative x.