6
$\begingroup$

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

$\endgroup$
4
  • 1
    $\begingroup$ You want to evaluate for $m$ from zero to $m$? Was that upper limit meant to be $M$? $\endgroup$ Jul 9, 2013 at 7:05
  • 1
    $\begingroup$ @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. $\endgroup$ Jul 9, 2013 at 11:17
  • $\begingroup$ @FelixGoldberg I thought about that, but the OP says $A$ and $B$ are positive-definite. $\endgroup$ Jul 9, 2013 at 12:56
  • $\begingroup$ @FelixGoldberg In mid/late 2018, I come across the Woodbury formula in a totally different setting, and that resulted in a flashback for this particular question several years back. Your Woodbury formula comment flew over my head initially, but now I am so glad that someone mentioned it originally. $\endgroup$
    – Yi Liu
    Aug 25, 2018 at 8:42

3 Answers 3

4
$\begingroup$

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.

$\endgroup$
3
$\begingroup$

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.

$\endgroup$
3
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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