MathOverflow is a question and answer site for professional mathematicians. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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

share|cite|improve this question
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
up vote 5 down vote accepted

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.

share|cite|improve this answer

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.

share|cite|improve this answer

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.

share|cite|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.