Fdtxjr

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)^-1x$ 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)^-1A^-1$

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

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.

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+1A^-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^-1A)^-1 = alpha^-1I - alpha^-2A + alpha^-3A - 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.)

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=fracdet(A)det(A+alpha I)A^-1+alpha I$$

Proof:

I'll use the formula of

$$ adj(A+alpha I) = adj(A) + \ + beginbmatrix
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) \
endbmatrix$$
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 = fracadj(A)det(A+alpha I) + alpha I$$
And using the $A^-1 det(A) = adj(A)$ again:
$$(A+alpha I)^-1=fracdet(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_11a_22...a_nn$$

if $A=[a_ij]_i,j in 1,2,...,n$.

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)^-1Q^T.$ Since $A$ is positive semidefinite, $(D+alpha I)^-1$ with $alpha>0$, exists even when $A^-1$ does not.