# 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 – Shogun Mar 26 at 3:42

## 3 Answers

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! – shahaf finder Mar 29 at 17:46
• What do you mean by GMM? – eepperly16 Mar 30 at 19:14
• The matrices are cov/precision matrices from a gaussian mixture model – shahaf finder Mar 30 at 19:23

## WRONG ANSWER (sorry, I did something wrong, please ignore!)

There is a way to write this up so that it contains $$det(A)$$, but not $$A$$:

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

Proof:

I'll use the formula of

$$adj(A+\alpha I) = adj(A) + \\ + \begin{bmatrix} \alpha \ det(A+\alpha I) & 0 & \dots & 0 \\ 0 & \alpha \ det(A+\alpha I) & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \alpha \ det(A+\alpha I) \\ \end{bmatrix}$$ Which comes from the definition of the Adjugate matrix: https://en.wikipedia.org/wiki/Adjugate_matrix

This is equivalent to: $$adj(A+\alpha I) = adj(A) + det(A+\alpha I)\ \alpha I$$ I'll use the equality of $$A^{-1} det(A) = adj(A)$$ and reorder the equation: $$det(A+\alpha I)(A+\alpha I)^{-1} = adj(A) + det(A + \alpha I)\ \alpha I$$ Diving both sides by $$det(A+\alpha I)$$ (if it's non-$$0$$): $$(A+\alpha I)^{-1} = \frac{adj(A)}{det(A+\alpha I)} + \alpha I$$ And using the $$A^{-1} det(A) = adj(A)$$ again: $$(A+\alpha I)^{-1}=\frac{det(A)}{det(A+\alpha I)}A^{-1}+\alpha I \quad \blacksquare$$

There might be an even simpler way of calculating $$(A+\alpha I)^{-1}$$, but this is the only thing that comes to my mind at the moment.

Edit: Also, $$det(A+\alpha I)$$ can be calculated from $$det(A)$$ and $$A$$'s main diagonal:

$$det(A+\alpha I) = det(A) + (a_{11} + \alpha)(a_{22} + \alpha)...(a_{nn} + \alpha) - a_{11}a_{22}...a_{nn}$$

if $$A=[a_{ij}]_{i,j \in \{1,2,...,n\}}$$.

• This sounds like a good option to consider, the only problem is the last part, because I need the diagonal values of A. Can I somehow get it without doing inverting $A^{-1}$? – shahaf finder Mar 25 at 18:53
• Well there is the equality of $det(A^{-1})=\frac{1}{det(A)}$, so $\frac{det(A)}{det(A+\alpha I)}=\frac{det((A+\alpha I)^{-1})}{det(A^{-1})}$. I can't think of anything else right now, but if I do, I'll comment about it. I also looked up the problem of the site and found this: math.stackexchange.com/questions/978051/… – Daniel P Mar 25 at 19:15
• $det(A)$ is not the problem, but $det(A+\alpha I)$ is. Ignoring this issue, I've done a quick test using python and it doesn't work, I'll keep on testing but there might be an error somewhere – shahaf finder Mar 25 at 19:20
• Where does the first formula come from? – shahaf finder Mar 25 at 19:51
• Thats ok, thanks for trying! – shahaf finder Mar 25 at 20:02

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}$. – shahaf finder Mar 29 at 13:49