Solving $(A^{-1} + D)^{-1} v$ with low rank Cholesky factors of $A$

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?

• 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? Sep 23 '15 at 17:21
• @GoHokies what kind of iterative solver do you have in mind? I'll be happy with any fast numerical solution. Sep 23 '15 at 18:11
• Related answers: scicomp.stackexchange.com/a/10631/1128 Sep 23 '15 at 20:36
• If $D$ is positive definite, I'd suggest you use CG on $(GG^T+D^{-1})$. Otherwise, SYMMLQ or MINRES. Sep 24 '15 at 8:58
• 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. Sep 24 '15 at 13:22

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.