Inverse of two matrix sums

Question

Trying to learn a basic implementation from the following post:

Inverse of the sum of matrices

and also from

https://www.jstor.org/stable/2690437?seq=1#page_scan_tab_contents

The post indicates that $$ (A+B)^{-1} = A^{-1} - \frac{1}{1+g}A^{-1}BA^{-1}. $$

where $g=\operatorname{trace}(BA^{-1})$.

The following is my R code implementation. Suppose I have two simple matrices

A <- structure(c(1, 0, 0, 1), .Dim = c(2L, 2L))
B <- structure(c(1.5, 0.5, 0.5, 0.5), .Dim = c(2L, 2L))

Then, the inverse of their sum is

> solve(A+B)
           [,1]       [,2]
[1,]  0.4285714 -0.1428571
[2,] -0.1428571  0.7142857

Now, I think the RHS of the example from above would be implemented as

> g <- sum(diag(B %*% A)) ### trace of matrix
> A - 1/(1+g) * A %*% B %*% A
           [,1]       [,2]
[1,]  0.5000000 -0.1666667
[2,] -0.1666667  0.8333333

I realize $A = I$, so this reduces even further. But, perhaps there is a portion of the lemma I'm misunderstanding given that the second implementation does not match the first.

In my real world problem, B is huge and taking its inverse inside an iterative algorithm would be costly. So, trying to find the cheapest implementation I can.

Any advice on how to implement this in an efficient way or correct my misunderstandings would be welcome.

Answer 1

The formula you quote is said to be valid if $rank(B)=1$.

And the issue is that it becomes invalid if $rank(B)>1$.

let us consider (as you do) the case $A=I$. The given formula would become

$$(I+B)^{-1}=I-aB \tag{1}$$

$$ \text{where} \ \ a:=\frac{1}{1+\operatorname{trace}(B)}\tag{2}$$

But formula (1) is false in general.

A counterexample among many: take $B=\begin{pmatrix}1&1\\0&1\end{pmatrix}$.

The LHS of (1) is: $\begin{pmatrix}0.5&-0.25\\0&0.5\end{pmatrix}$

while its RHS is $\begin{pmatrix}0.8&-0.2\\0& \ \ 0.8\end{pmatrix}.$