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

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.

-

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.

-

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.

-