Fix $n \times n$ real symmetric positive definite matrices $A$ and $B$. Fix vectors $x$ and $y$ in $\mathbf{R}^n$. I want to compute the following bilinear products quickly: $\{x^T (A+mB)^{-1} y\}\_{m=0}^M$.


A naive but practical method would involve inverting each matrix individually, thus requiring $O(Mn^3)$.


Are there practical improvements for $M > n > \log(M)$ ? (e.g. Hopefully $O(n^3 + n^2M\log(n))$ or something).



Algorithm to quickly compute the individual inverses of a linear sequence of matrices
Here is a faster way (perhaps not the nicest from a numerical stability point of view).


\begin{eqnarray\*}
x^T(A+mB)^{-1}y &=& x^T(A^{1/2}(I+mA^{-1/2}BA^{-1/2})A^{1/2})^{-1}y\\
&=& x^TA^{-1/2}(I+mZ)^{-1}A^{-1/2}y\\
&=& p^T(I+mZ)^{-1}q.
\end{eqnarray\*}
Now let $Z=U\Lambda U^T$, so that we need to compute 
\begin{equation\*}
p^TU(I+m\Lambda)^{-1}U^Tq = \tilde{p}^T(I+m\Lambda)^{-1}\tilde{q},\qquad m=0,1,...,M.
\end{equation\*}
It takes $O(n^3)$ to compute $Z$, the eigendecomposition of $Z$, as well as to compute vectors $p$ and $q$. The overall cost is thus $O(n^3 + Mn)$, as desired.


