5
$\begingroup$

I have a large matrix $A \in \mathcal{R}^{N\times N}$ which is supposedly positive-definite, but numerically low rank. Instead of $A$, I have its incomplete Cholesky factor $G$, such that $A \simeq GG^\top$. $A$ is very large, so I can't keep it in memory, but $G$ is small enough. I need to numerically solve:

$$ (A^{-1} + D)^{-1} v $$

where $D$ is a diagonal matrix with positive entries, and $v$ is a vector. I want to do so without forming an $N \times N$ matrix. Is this possible without further approximation?

(Context: this is the Newton-step in an optimization for variational Bayes under GP. $A$ is a Gram matrix, which I can evaluate any value if needed.)

Using the matrix inversion lemma, we could rewrite it as,

$$ (A - A (D^{-1} + A)^{-1} A)m \simeq G(G^\top m) - GG^\top (D^{-1} + A)^{-1} G(G^\top m)$$

But, the $(D^{-1} + A)^{-1}$ term is still troublesome. (Perhaps a low-rank approximation of it could help?)

EDIT: If only I had the incomplete Cholesky of $A^{-1} = L L^\top$, things would have been much easier. In this case, matrix inversion lemma gives,

$$ (D^{-1} - D^{-1} L (I - L^\top D L)^{-1} L^\top D^{-1}) m $$

which can be all computed in the low-rank dimensions. Now, note that the pseudo-inverse of $G$ gives $L$. Hence, an econ-SVD on $G$ would solve it. Is there a better way?

Improve this question
$\endgroup$
7
  • 3
    $\begingroup$ Would an iterative solver count as "further approximation"? Even if the relative error of the linear solver is kept smaller than the "truncation errors" incurred in the incomplete Cholesky decomposition? $\endgroup$
    – GoHokies
    Sep 23, 2015 at 17:21
  • $\begingroup$ @GoHokies what kind of iterative solver do you have in mind? I'll be happy with any fast numerical solution. $\endgroup$
    – Memming
    Sep 23, 2015 at 18:11
  • $\begingroup$ Related answers: scicomp.stackexchange.com/a/10631/1128 $\endgroup$
    – Memming
    Sep 23, 2015 at 20:36
  • 1
    $\begingroup$ If $D$ is positive definite, I'd suggest you use CG on $(GG^T+D^{-1})$. Otherwise, SYMMLQ or MINRES. $\endgroup$
    – GoHokies
    Sep 24, 2015 at 8:58
  • 2
    $\begingroup$ Try Conjugate Gradients (CG). With CG (or other Lanczos-type iterative solvers) you do not need to form $GG^T + D^{-1}$ explicitly, just to calculate its action on a vector: $(GG^T + D^{-1}) y = G(G^T y) + D^{-1}y$. Total (asymptotic) cost of such a vector product: ${\cal O}(N^2)$ operations (the real cost is even lower if your Cholesky factor is sparse). One further observation: your linear solve is the inner loop of a Newton iteration, so it may pay off (computationally) to do less exact solves at the beginning, and tighten the tol-factor as you approach the optimal (Newton) solution. $\endgroup$
    – GoHokies
    Sep 24, 2015 at 13:22

1 Answer 1

Reset to default
2
$\begingroup$

Y.Z. found a simpler solution without an explicit pseudo-inverse.

Using the matrix inversion lemma twice, he got: $$(A^{-1} + D)^{-1}v = GG^\top v - GG^\top D GG^\top v + GG^\top D G (I + G^\top D G)^{-1} G^\top D GG^\top v$$ which can be evaluated without forming the large $N \times N$ matrix.

Improve this answer
$\endgroup$

Your Answer

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

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