2

I need to calculate $(A+B)^{-1}$, where $A$ and $B$ are two square, very sparse and very large. $A$ is block diagonal, real symmetric and positive definite, and I have access to $A^{-1}$ (which in this case is also sparse, and block diagonal). $B$ is diagonal and real positive. In my application, I need to calculate the inverse of the sum of these two matrices where the inverse of the non-diagonal one (e.g. $A^{-1}$) is updated frequently, and readily available.

Since $B$ is full rank, the Woodbury lemma is of no use here (well, it is, but it's too slow). Other methods described in this nice question are of no use in my case as the spectral radius of $A^{-1}B$ is much larger than one. Methods based on a diagonalisation assume that it is the diagonal matrix that is being updated frequently, whereas that's not my case (i.e., diagonalising $A$ is expensive, and I'd have to do that very often).

I'm quite happy to live with an approximate solution.

Jose
  • 131

2 Answers2

4

Let $A_1,A_2,\cdots,A_q$ be the diagonal blocks of $A$, and $a_{1,1},a_{1,2},\cdots,a_{1,n_1},a_{2,1},a_{2,2},\cdots,a_{2,n_2},\cdots,a_{q,1},a_{q,2},\cdots,a_{1,n_q}$ the diagonal elements of $B$, then the inverse of the sum would simply be a diagonal block matrix with blocks: ${(A_i+diag(a_{i,1},\cdots,a_{i,n_i}))}^{-1}$ for $i\in(1,2,\cdots,q)$.

So the problem is reduced to finding the inverse of the sum of a matrix and a diagonal matrix. Fortunately, $A$ is symmetric positive definite, so each $A_i$ diagonalizable, hence, we can write it as follows:

$$ A_i=P_iD_i{P_i}^{T}=P_i\begin{bmatrix} \lambda_{i,1} & 0 & 0 & \cdots & 0 & 0 \\ 0 & \ddots & \ddots & \ddots & & 0 \\\ 0 & \ddots & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & 0 \\ 0 & & \ddots & \ddots & \ddots & 0 \\ 0 & 0 & \cdots & 0 & 0 & \lambda_{i,n_{i}} \end{bmatrix}{P_i}^{T} $$

Where $\lambda_{i,1},\cdots,\lambda_{i,n_{i}}$ are the eigenvalues of $A_i$, hence:

$$ A_i+diag(a_{i,1},\cdots,a_{i,n_i})=P_i\begin{bmatrix} \lambda_{i,1}+a_{n_i} & 0 & 0 & \cdots & 0 & 0 \\ 0 & \ddots & \ddots & \ddots & & 0 \\\ 0 & \ddots & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & 0 \\ 0 & & \ddots & \ddots & \ddots & 0 \\ 0 & 0 & \cdots & 0 & 0 & \lambda_{i,n_{i}}+a_{i,n_i} \end{bmatrix}{P_i}^{T} $$

Since $A_i$ and $B$ are symmetric positive definite, then $\lambda_{i,j}+a_{i,j}\neq 0$, so :

$$ {(A_i+diag(a_{i,1},\cdots,a_{i,n_i}))}^{-1}=P_i\begin{bmatrix} \frac{1}{\lambda_{i,1}+a_{i,1}} & 0 & 0 & \cdots & 0 & 0 \\ 0 & \ddots & \ddots & \ddots & & 0 \\\ 0 & \ddots & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & 0 \\ 0 & & \ddots & \ddots & \ddots & 0 \\ 0 & 0 & \cdots & 0 & 0 & \frac{1}{\lambda_{i,n_{i}}+a_{i,n_i}} \end{bmatrix}{P_i}^{T}=P_iD_i{P_i}^{T} $$

Hence the inverse of ${(A+B)}^{-1}$ is a block diagonal matrix with diagonal elements being the matrices above.

You can rewrite that as

$$ {(A+B)}^{-1}=\begin{bmatrix} P_1 & 0 & 0 & \cdots & 0 & 0 \\ 0 & \ddots & \ddots & \ddots & & 0 \\\ 0 & \ddots & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & 0 \\ 0 & & \ddots & \ddots & \ddots & 0 \\ 0 & 0 & \cdots & 0 & 0 & P_q \end{bmatrix}\times\begin{bmatrix} D_1 & 0 & 0 & \cdots & 0 & 0 \\ 0 & \ddots & \ddots & \ddots & & 0 \\\ 0 & \ddots & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & 0 \\ 0 & & \ddots & \ddots & \ddots & 0 \\ 0 & 0 & \cdots & 0 & 0 & D_q \end{bmatrix}\times\begin{bmatrix} {P_1}^{T} & 0 & 0 & \cdots & 0 & 0 \\ 0 & \ddots & \ddots & \ddots & & 0 \\\ 0 & \ddots & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & 0 \\ 0 & & \ddots & \ddots & \ddots & 0 \\ 0 & 0 & \cdots & 0 & 0 & {P_q}^{T} \end{bmatrix} $$

Oussama Boussif
  • 3,085
0

Unless I am seriously mistaken, the above answer is incorrect. In my opinion, the problem lies with $$ A_i+diag(a_{i,1},\cdots,a_{i,n_i})=P_i\begin{bmatrix} \lambda_{i,1}+a_{n_i} & 0 & 0 & \cdots & 0 & 0 \\ 0 & \ddots & \ddots & \ddots & & 0 \\\ 0 & \ddots & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & 0 \\ 0 & & \ddots & \ddots & \ddots & 0 \\ 0 & 0 & \cdots & 0 & 0 & \lambda_{i,n_{i}}+a_{i,n_i} \end{bmatrix}{P_i}^{T} $$

Let's consider $$ P=\frac{1}{\sqrt{5}}\left[\begin{matrix} 2 & -1 \\ 1 & 2\end{matrix}\right] $$ $$ D=\left[\begin{matrix} 2 & 0 \\ 0 & 1\end{matrix}\right] $$ $$ B= \left[\begin{matrix} 2 & 0 \\ 0 & 1\end{matrix}\right] $$

Quick calculations show $$ A + B = PAP^T + B = \left[\begin{matrix} 3.8 & 0.4 \\ 0.4 & 2.2\end{matrix}\right], $$ but $$ P ( D + B) P' = \left[\begin{matrix} 3.6 & 0.8 \\ 0.8 & 2.4\end{matrix}\right]. $$

The problem lies in the fact that adding an element to the diagonal of an un-transformed matrix A is not the same as adding it to the diagonal of the diagonalization matrix, i.e. $$P A P^T + B \neq P (A + B) P^T$$ in general.

Or am I missing something?

PS: Unfortunately, I do not have enough reputation to comment. Sorry for the second answer, please feel free to add this to the original.

theo
  • 21
  • 3