Given real-valued vectors $\mathbf{a},\mathbf{u}$ in $\mathbb{R}^d$, is there a nice expression the inverse Laplace transform of $f(y)$ below?
$$f(y)=\frac{\left(\sum_i \frac{u_i}{y-a_i}\right)^2}{1-\sum_i u_i \frac{u_i}{y-a_i}}$$
If nice formula for $g(t)=\mathscr{L}^{-1}[f]$ is not possible, I'm interested in a nice upper bound on $g(t)$ which works for my application. I have $a=-2u(1-u)$ and $u$ represents a discrete probability distribution over $d$ outcomes (entries of $u$ are positive, adding up to 1)
Motivation:
$f(y)$ is the Laplace transform of $g(t)=\langle \exp(t(A+B))-\exp(tA),\mathbf{1}\rangle$ where $A$ is diagonal matrix with $\mathbf{a}$ on diagonal and $B=\mathbf{u}\mathbf{u}^T$ is rank-1. Inverting this Laplace transform gives a method to compute rank-1 correction to matrix exponential and matrix power.
Mathematica seems to discover simple formulas for $g$ for specific values of $\mathbf{a}$,$\mathbf{u}$, I'm looking for a more generic expression. In my application, $\mathbf{a}=-2 \mathbf{u} (1 - \mathbf{u})$ while entries of $\mathbf{u}$ are strictly positive and adding up to 1.
For instance, it tells me that
$$L^{-1}\text{[f](t)=}-1. e^{-0.495868 t}+0.112202 e^{-0.423048 t}-1. e^{-0.396694 t}+0.0556054 e^{-0.317963 t}-1. e^{-0.297521 t}+2.83219 e^{-0.0441126 t}$$ when
$$\mathbf{a}=\left( \begin{array}{c} -\frac{60}{121} \\ -\frac{48}{121} \\ -\frac{36}{121} \\ \end{array} \right)$$
$$\mathbf{u}=\left( \begin{array}{c} \frac{6}{11} \\ \frac{3}{11} \\ \frac{2}{11} \\ \end{array} \right)$$
If we do the partial fraction decomposition of your Laplace transform we get:
\begin{eqnarray} f[y]&=& - \sum_{j=1}^d \frac{1}{y -a_j} + \frac{\sum\limits_{\xi=0}^{d-1} {\mathfrak B}_\xi^{(d)} y^\xi} {\sum\limits_{\xi=0}^d {\mathfrak A}_\xi^{(d)} y^\xi} \\ &=& - \sum_{j=1}^d \frac{1}{y -a_j} + \sum\limits_{\xi=0}^{d-1} \frac{{\mathfrak B}_\xi^{(d)}}{{\mathfrak A}^{(d)}_d} \frac{y^\xi}{ \prod\limits_{\eta=1}^d (y - \zeta_\eta)} \\ &=& - \sum_{j=1}^d \frac{1}{y -a_j} + \sum\limits_{\xi=0}^{d-1} \frac{{\mathfrak B}_\xi^{(d)}}{{\mathfrak A}^{(d)}_d} \sum\limits_{\eta=1}^d \frac{(\zeta_\eta)^\xi}{y - \zeta_\eta} \cdot \frac{1}{\prod\limits_{\lambda=1,\lambda\neq \eta}^d (\zeta_\eta - \zeta_\lambda)} \end{eqnarray}
where $(\zeta_\eta)_{\eta=1}^d $ are roots of the polynomial in the denominator in the first equation on the top.
Now, I took ${\bf a} = -2 {\bf u}(1- {\bf u})$ and, as you wrote, ${\bf u}$ being strictly positive and summing up to one and the roots seem to be distinct and always real, at least for $d=2,\cdots,5$, as seen from the code snippet below:
Then inverting the expression in the last equation is straightforward because ${\mathfrak L}^{-1}_y [ 1/(y+a)] (t) = e^{-a t}$.
Update:
The polynomial in the denominator in the right hand side reads:
\begin{equation} \sum\limits_{\xi=0}^d {\mathfrak A}_\xi^{(d)} y^\xi = \prod\limits_{j=1}^d (y-a_j) - \sum\limits_{i=1}^d u_i^2 \prod\limits_{j=1,j\neq i}^d (y-a_j) \end{equation} and ${\mathfrak A}_d^{(d)} = 1$ whereas the coefficients $\left({\mathfrak B}_\xi^{(d)}\right)_{\xi=0}^{d-1}$ read:
\begin{eqnarray} \left( \begin{array}{c} B_0\to -a_1-a_2-\left(u_1-u_2\right){}^2 \\ B_1\to 2 \\ \end{array} \right) \end{eqnarray} for $d=2$.
\begin{eqnarray} \left( \begin{array}{c} B_0\to a_3 \left(u_1-u_2\right){}^2+a_2 \left(a_3+\left(u_1-u_3\right){}^2\right)+a_1 \left(a_2+a_3+\left(u_2-u_3\right){}^2\right) \\ B_1\to -2 \left(a_1+a_2+a_3+u_1^2-u_2 u_1-u_3 u_1+u_2^2+u_3^2-u_2 u_3\right) \\ B_2\to 3 \\ \end{array} \right) \end{eqnarray} for $d=3$.
\begin{eqnarray} \left( \begin{array}{c} B_0\to -a_3 a_4 \left(u_1-u_2\right){}^2-a_2 \left(a_4 \left(u_1-u_3\right){}^2+a_3 \left(a_4+\left(u_1-u_4\right){}^2\right)\right)-a_1 \left(a_4 \left(u_2-u_3\right){}^2+a_3 \left(a_4+\left(u_2-u_4\right){}^2\right)+a_2 \left(a_3+a_4+\left(u_3-u_4\right){}^2\right)\right) \\ B_1\to 2 \left(a_3 u_1^2+a_4 u_1^2-a_3 u_2 u_1-a_4 u_2 u_1-a_4 u_3 u_1-a_3 u_4 u_1+a_3 u_2^2+a_4 u_2^2+a_4 u_3^2+a_3 u_4^2-a_4 u_2 u_3-a_3 u_2 u_4+a_2 \left(a_3+a_4+u_1^2-u_3 u_1-u_4 u_1+u_3^2+u_4^2-u_3 u_4\right)+a_1 \left(a_2+a_3+a_4+u_2^2-u_3 u_2-u_4 u_2+u_3^2+u_4^2-u_3 u_4\right)+a_3 a_4\right) \\ B_2\to -3 a_1-3 a_2-3 a_3-3 a_4-3 u_1^2+2 u_2 u_1+2 u_3 u_1+2 u_4 u_1-3 u_2^2-3 u_3^2-3 u_4^2+2 u_2 u_3+2 u_2 u_4+2 u_3 u_4 \\ B_3\to 4 \\ \end{array} \right) \end{eqnarray} for $d=4$.
For higher values of $d \ge 5$ those coefficients can be generated using the code snippet below:
Question:
Why is that inverse Laplace transform equal to that quantity that involves matrix exponentials? Can you explain it?