Inverse of the sum of a symmetric and diagonal matrices

This is more a long comment or a remark than an answer, but maybe it can help, at least in certain cases.

Assume $B$ is invertible, and so is $A+\delta B$.

Notice that $$A+\delta B = \left(\frac{1}{\delta}AB^{-1}+1\right)\delta B.$$ Therefore $$(A+\delta B)^{-1} = \frac{1}{\delta}B^{-1}\left(\frac{1}{\delta}AB^{-1}+1\right)^{-1}$$ Now assume that $\left|\det\left(\frac{1}{\delta}AB^{-1}\right)\right|<1$. Then we have $$\left(\frac{1}{\delta}AB^{-1}+1\right)^{-1} = \sum_{n\ge0}(-\delta)^{-n}(AB^{-1})^n.$$ In conclusion, we have $$(A+\delta B)^{-1} = -\sum_{n\ge0}(-\delta)^{-(n+1)}B^{-1}(AB^{-1})^n.$$

Assume that $\delta$ varies in a set with $p$ elements and $A,B\in M_n(\mathbb{R})$. Then the complexity of the trivial calculation of the $(A+\delta B)^{-1}$ is $O(pn^3)$. Implicitly you write that we can do better when $B=I_n$; do you want to bet ?

Point of view 1. We assume that $A,A+\delta B$ are invertible and that $A^{-1}B$ is diagonalizable: $A^{-1}B=PDP^{-1}$. Let $(\lambda_i)_i=spectrum(A^{-1}B)$. Then $(A+\delta B)^{-1}=A^{-1}-\delta A^{-1}B(I+\delta A^{-1}B)^{-1}A^{-1}=A^{-1}- A^{-1}BPdiag(\mu_1,\cdots,\mu_n)P^{-1}A^{-1}$, where $\mu_i=\delta/(1+\delta\lambda_i)$. After precalculation of $P,D,U=A^{-1}BP,V=P^{-1}A^{-1}$, the required complexity is that of calculation of $Udiag(\mu_1,\cdots,\mu_n)V$ (a tensor product). If we stack the result matrix into a column, then there is a fixed $W\in M_{n^2 ,n}$ s.t. the calculation reduces to $W[\mu_1,\cdots,\mu_n]^T$, the complexity of which, is $O(n^3)$.

If $B=I, spectrum(A)=(\alpha_i)_i$ and $A$ is diagonalizable: $A=QD_1Q^{-1}$, then $(A+\delta I)^{-1}=Qdiag(\nu_i)Q^{-1}$, where $\nu_i=1/(\alpha_i+\delta)$. Again, the complexity is $O(n^3)$.

Point of view 2. $(A+\delta B)^{-1}$ is a rational fraction of $\delta$ in the form $\dfrac{1}{p_n(\delta)}Z$ where $z_{i,j}(\delta)$ is a polynomial of degree $\geq n-2$ and $p_n$ is a polynomial of degree $n$ (when $B=I$, the form is not simpler!). The required complexity is the complexity of the evaluation of $n^2$ polynomials of degree $\geq n-2$. With the Horner method, the complexity is $n^2$ times $(n$ mult, $n$ add). Using a precalculation and the Knuth/Eve method, we can do better: $n^2$ times $(n/2$ mult, $n$ add); thus, the complexity is almost halved (in fact, it is not really true because the gap between the complexity of multiplication and addition has decreased lately).

We can also write the result in the form (with precalculation of course) $(A+\delta B)^{-1}=\dfrac{1}{p_n(\delta)}(U_0+\delta U_1+\cdots \delta^{n-2}U_{n-2}+\delta^{n-1}\Delta)$ where $U_i,\Delta\in M_n$ and $\Delta$ is diagonal; again, the associated complexity is $O(n^3)$.

Conclusion: do not expect miracles.

EDIT. I just read your incomprehensible edit. That is clear: even, after precalculations, you cannot (eventually) diagonalize $A+\delta B$ with complexity $o(n^3)$ (except if $B$ is a homothety).