Assume I want to model something that has happened as a diffusion process and I have two points in time $t=0, t=N$ and I want to model what has happened to a function $f_0(x) \to f_N(x)$ in a way so that the "same thing" has occured $N$ times. In our case I imagine a convolution: $$f_{n+1}(x) = \int_{-\infty}^{\infty} f_{n}(\xi)g(x-\xi)d\xi$$
We know $f_0(x),f_N(x)$ and we want to calculate $g(x)$. I also if possible want the capability to push $g(x)$ to be as "lazy" as possible in some suitable sense - at least in case many different solutions would be feasible (as often is the case with inverse problems).
How can I formulate this mathematically? And how to solve for it?
Own work my own ideas so far are centered around discretizing $f_0(x)$, $f_N(x)$ on a grid, then do the same for unknown $f_1(x),f_2(x),\cdots ,f_{N-1}(x)$. Then start with some "guess" of how these $f_k(x)$ behave. Then to stuff their discretizations into a looong vector, call it $\bf f = [f_0,f_1,\cdots,f_N]^T$, then call the matrix representing $(\cdot)*g(x)$ for $\bf G$ and pose the following equation system:
$$\left[\begin{array}{cccccc}\bf G&\bf -I&\bf 0&\cdots&\bf 0&\bf 0\\\bf 0&\bf G&\bf -I&\bf 0&\cdots&\bf 0\\\bf 0&\ddots&\bf G&\bf -I&\ddots&\bf 0\\\bf 0&\ddots&\ddots&\ddots&\ddots&\ddots\\\bf 0&\bf 0&\bf 0&\bf 0&\bf G&\bf -I\end{array}\right]\left[\begin{array}{c}\bf f_0\\\bf f_1\\\bf f_2\\\vdots\\\bf f_{N-1}\\\bf f_N\end{array}\right] = \left[\begin{array}{c}\bf 0\\\bf 0\\\bf 0\\\vdots\\\bf 0\\\bf 0\end{array}\right]$$
Then to alternatingly solve this for $\bf G$ (with suitable regularization) and for $\bf f$ (with boundary conditions and regularization).
If we solve the equation system above for $\bf G$ and $\bf f$, then $${\bf G}^N {\bf f_0 = f_N}$$, right?
Yes, this is doable, even quite impressive results can be had by simplifying to $$\min_{\bf G}\left\{\sum_k \|{\bf [G,-I,0,\cdots,0]C}^k {\bf f}\|_F^2\right\}$$ for the cyclic permutation operator $\bf C$, encouraging $\bf G$ to be Toeplitz with an extra cost-term and of course other regularizing terms like punishing "filter taps" too far off the diagonal. This we can achieve by adding a cost on the convolution of $\left[\begin{array}{rr}1&0\\0&-1\end{array}\right]$ on $\bf G$ as if it were a discretized 2D function fulfilling a differential equation. We use only one naiive guess of a sequence of $$f_k(t) = \frac k N f_N(t) + \left(1-\frac k N\right)f_0(t)$$ and can get results like this, even without any iteration, solving for better $f_k(t)$ candidate functions.