Efficient diagonal update of matrix inverse

I am computing $(kI + A)^{-1}$ in an iterative algorithm where $k$ changes in each iteration. $I$ is an $n$-by-$n$ identity matrix, $A$ is an $n$-by-$n$ precomputed symmetric positive-definite matrix. Since $A$ is precomputed I may invert, factor, decompose, or do anything to $A$ before the algorithm starts. $k$ will converge (not monotonically) to the sought output.

Now, my question is if there is an efficient way to compute the inverse that does not involve computing the inverse of a full $n$-by-$n$ matrix?

• Turns out this was very simple. I'll just write it here if someone else has need for it: $(kI+A)^{-1}=(kI+PDP^{-1})^{-1}=(P(D + kI)P^{-1})^{-1}=P(D + kI)^{-1}P^{-1}$, where $A=PDP^{-1}$ is the eigenvalue decomposition. And the inverse of a diagonal matrix is quickly computed as the matrix with diagonal elements the reciprocal of the diagonal elements of the original matrix. This is roughly 50 times faster than the original solution for my data on my computer. Feb 17, 2014 at 14:22
• I have a similar question, math.stackexchange.com/q/2667566/10063. But my update of the diagonal has distinct entries along the diagonal, so your method doesn't work.
– a06e
Feb 26, 2018 at 15:16

Turns out this was very simple. I'll just write it here if someone else has need for it: $$(kI+A)^{-1}=(kI+PDP^{-1})^{-1}=(P(D + kI)P^{-1})^{-1}=P(D + kI)^{-1}P^{-1}$$ where $A=PDP^{-1}$ is the eigenvalue decomposition. And the inverse of a diagonal matrix is quickly computed as the matrix with diagonal elements the reciprocal of the diagonal elements of the original matrix. This is roughly $50$ times faster than the original solution for my data on my computer.

-- Tommy L

EDIT. 1. The Tommy L 's method is not better than the naive method.

Indeed, the complexity of the calculation of $$(kI+A)^{-1}$$ is $$\approx n^3$$ blocks (addition-multiplication).

About the complexity of $$P(D+kI)^{-1}P^{-1}=QP^{-1}$$ (when we know $$P,D,P^{-1}$$); the complexity of the calculation of $$Q$$ is $$O(n^2)$$ and the one of $$QP^{-1}$$ is $$\approx n^3$$ blocks as above. Note that one works with a fixed number of digits.

1. When the real $$k$$ can take many (more than $$n$$) values, one idea is to do the following

The problem is equivalent to calculate the resolvent of $$A$$: $$R(x)=(xI-A)^{-1}=\dfrac{Adjoint(xI-A)}{\det(xI-A)}=$$

$$\dfrac{P_0+\cdots+P_{n-1}x^{n-1}}{a_0+\cdots+a_{n-1}x^{n-1}+x^n}$$.

STEP 1. Using the Leverrier iteration, we can calculate the matrices $$P_i$$ and the scalars $$a_j$$

$$P_{n-1}:=I:a_{n-1}:=-Trace(A):$$

for $$k$$ from $$n-2$$ by $$-1$$ to $$0$$ do

$$P_k:=P_{k+1}A+a_{k+1}I:a_k:=-\dfrac{1}{n-k}Trace(P_kA):$$

od:

During this step, we must calculate the exact values of the $$P_i,a_j$$ (assuming that $$A$$ is exactly known), which requires a very large number of digits. Then the complexity of the calculation is $$O(n^4)$$ (at least) but it's done only one time.

STEP 2. We put $$x:=-k_1,-k_2,\cdots$$. Unfortunately, the time of calculation -with a fixed number of significative digits- of $$R(-k)$$, is larger than the time of calculation of $$(-kI-A)^{-1}$$ !!!