# Inverse of SPD matrix + identity

I didn't know exactly where to post this question. I feel like it falls between Computer Science and Mathematics, so I'm sorry if it doesn't fit here.

I need to calculate $$(A+\alpha I)^{-1}$$, given $$\alpha>0$$ and $$A^{-1}$$ which is SPD, and with a known sparsity pattern. (If it helps in any way, I need it for calculating Mahalanobis distance)

I'm dealing with an high dimension and sparse $$A^{-1}$$ so I would also like to avoid calculating $$A$$ (or any other inverse) using the inverse operation.

I tried looking into Woodbury Matrix Identity, but I can't find a way to use it in my case.
Is there any closed form solution or iterative method that I can use?
Is the fact that I need only to calculate $$x^T(A+\alpha I)^{-1}x$$ can help in any way?

update:
I found an interesting way to avoid calculating $$A$$ out of $$A^{-1}$$ for this:

$$(A+\alpha I)^{-1} = (A+\alpha AA^{-1})^{-1} = (A(I+\alpha A^{-1}))^{-1} = (I+\alpha A^{-1})^{-1}A^{-1}$$

So now when calculating the Mahalanobis distance I need:
$$x^T(I+\alpha A^{-1})^{-1}A^{-1}x$$

Now I only need to do one inverse operation.
$$A^{-1}$$ is somewhat of a k-diagonal matrix.
So maybe now I'll find a way to calculate what i need more efficiently.

• Look up krylov methods
– user3417
Mar 26, 2019 at 3:42

As you suggest, the Woodbury Matrix Identity is of no use since the perturbation of $$\alpha I$$ is not low rank. In addition, in general, computing a full eigendecomposition $$A=QDQ^T$$ will be much slower than just using sparse Gaussian elimination on $$A+\alpha I$$, so this won't be helpful either in practice.

As a commenter suggested, this may be a good time for preconditioned conjugate gradient, a Krylov subspace method. For $$\alpha$$ small, we have $$A + \alpha I \approx A$$ so $$A^{-1}(A+\alpha I) \approx I$$. Thus, $$A^{-1}$$, which you already know, should provide a good preconditioner. If a better preconditioner is needed and $$\alpha \|A^{-1}\| < 1$$ for an appropriate matrix norm, then we have the Neumann series

$$(A+\alpha I)^{-1} = A^{-1}(I+\alpha A^{-1})^{-1} = A^{-1} - \alpha A^{-2} + \alpha^2 A^{-3}-\cdots.$$

You can truncate this infinite series up to the $$A^{-k}$$ power and evaluate the preconditioner $$(A^{-1} - \alpha A^{-2} + \alpha^2 A^{-3} -\cdots + (-1)^{k+1}A^{-k})x$$ using $$k$$ matrix-vector products with $$A^{-1}$$ using Horner's method.

If instead $$\alpha$$ is very large, specifically $$\|A\|/\alpha < 1$$, then we can instead use the Neuman series

$$(A+\alpha I)^{-1} = \alpha^{-1}(I+\alpha^{-1}A)^{-1} = \alpha^{-1}I - \alpha^{-2}A + \alpha^{-3}A - \cdots$$

as a preconditioner. If $$\alpha$$ is somewhere in between, then neither of these Neumann series ideas will work, and you might want to investigate another preconditioner if you intend to use Conjugate Gradient. (Algebraic multigrid may work well, for instance.)

• Thanks for the answer. All my $A^{-1}$ come from a special case GMM, so I don't know if I can bound $||A||/\alpha$. But I'll check it! Mar 29, 2019 at 17:46
• What do you mean by GMM? Mar 30, 2019 at 19:14
• The matrices are cov/precision matrices from a gaussian mixture model Mar 30, 2019 at 19:23

For every symmetric real matrix $$A$$ there exists a real orthogonal matrix $$Q$$ such that $$Q^TAQ$$ is a diagonal matrix. Every symmetric matrix is thus, up to choice of an orthonormal basis, a diagonal matrix. If you can find it, then $$A=QDQ^T$$ and your expression becomes $$Q(D+\alpha I)^{-1}Q^T.$$ Since $$A$$ is positive semidefinite, $$(D+\alpha I)^{-1}$$ with $$\alpha>0$$, exists even when $$A^{-1}$$ does not.

• Thanks for the answer, A is positive definite, so $A^{-1}$ exists (and in fact it's whats given). The complexity of eigendecomposition is also $O(n^3)$ which is like inverting (if I remember correctly). I'll update the post with something I found to somewhat reduce the calculation of $A$ out of $A^{-1}$. Mar 29, 2019 at 13:49