# 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).

-
You want to evaluate for $m$ from zero to $m$? Was that upper limit meant to be $M$? –  Gerry Myerson Jul 9 '13 at 7:05
@FedericoPoloni's solution is very nice. If $B$ happens to have nice structure (like low-rank) then even better solutions may be possible, using the Woodbury formula. –  Felix Goldberg Jul 9 '13 at 11:17
@FelixGoldberg I thought about that, but the OP says $A$ and $B$ are positive-definite. –  Federico Poloni Jul 9 '13 at 12:56

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.