Algorithm to quickly compute the individual inverses of a linear sequence of matrices 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).
 A: 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.
A: The quantity you wish to compute is a rational function in $m$ (if you consider it as a formal unknown), of degree bounded by $n$. You can get its coefficients via rational interpolation by evaluating it for $2n$ different values of $m$ (plus or minus 1, I didn't bother getting the exact number here). Then you just evaluate the obtained rational function for any other value of $M$.
This should give you an algorithm with complexity $O(n^4+nM)$ and somewhat dubious stability.
A: The following ideas are somewhat similar to that in Suvrit's answer.
Idea 1. Instead of first computing the square root of $A$, just solve the generalized eigenvalue problem $Ax = \lambda B x$ directly. This is especially tractable since according to the OP both $A$ and $B$ are symmetric and positive definite. You get a matrix $U$ such that $U^T B U = I$, $U^T A U = \mathrm{diag} (\lambda_i)$, and
\begin{align}
  U^T (A + mB) U &= \mathrm{diag} (\lambda_i + m) , \\
  U^{-1} (A + mB)^{-1} (U^{-1})^T &= \mathrm{diag} (\frac{1}{\lambda_i + m}) .
\end{align}
The formula that you want is then
$$
  x^T (A + mB)^{-1} y = x^T U \mathrm{diag} (\frac{1}{\lambda_i + m}) U^{T} y .
$$
The cost of solving the generalized eigenvalue problem is still $O(n^3)$ and the cost of each of the above multiplications once $U$ is known is $O(n^2)$. Thus, the overall cost is $O(n^3 + Mn^2)$.
Idea 2. Now, an eigenvalue problem requires iteration. If you'd prefer an algorithm with a definite number of steps, you could instead simultaneously tridiagonalize $A$ and $B$, that is, find a matrix $U$ such that $U^T A U = T_A$ and $U^T B U = T_B$, where both $T_A$ and $T_B$ are tridiagonal. Tridiagonalization is in any case a common preconditioning step in eigenvalue problem algorithms. Similar to above, the formula that you want is then
$$
  x^T (A + mB)^{-1} y = x^T U (T_A + m T_B)^{-1} U^{T} y .
$$
Since $T_A + m T_B$ is tridiagonal, you can compute its inverse in linear time, though slower by a constant factor than a diagonal matrix.
The cost of a tridiagonalization is once again $O(n^3)$ and the cost of each of the above inversions + multiplications is still $O(n^2)$ (dense matrix-vector multiplication dominates tridiagonal inversion). Thus, the overall cost is again $O(n^3 + Mn^2)$, though without requiring iteration.
