Inverse of constant matrix plus diagonal matrix

Question

Is there an efficient way to calculate the inverse of an $N \times N$ diagonal matrix plus a constant matrix? I am looking at $N$ of around $40,000$.

$$\left[\begin{array}{cccc} a & b & \cdots & b\\ b & a & & \vdots\\ \vdots & & \ddots & b\\ b & \cdots & b & a \end{array}\right]^{-1} = \,\,?$$

Putting this in to mathematica, for $N \in \{2, 3, 4\}$, the result is:

$$\left[ \begin{array}{cc} a & b \\ b & a \\ \end{array} \right]^{-1} = \left[ \begin{array}{cc} \frac{a}{a^2-b^2} & -\frac{b}{a^2-b^2} \\ -\frac{b}{a^2-b^2} & \frac{a}{a^2-b^2} \\ \end{array} \right]$$

$$\left[ \begin{array}{ccc} a & b & b \\ b & a & b \\ b & b & a \\ \end{array} \right]^{-1} = \left[ \begin{array}{ccc} \frac{a^2-b^2}{a^3-3 a b^2+2 b^3} & \frac{-a b+b^2}{a^3-3 a b^2+2 b^3} & \frac{-a b+b^2}{a^3-3 a b^2+2 b^3} \\ \frac{-a b+b^2}{a^3-3 a b^2+2 b^3} & \frac{a^2-b^2}{a^3-3 a b^2+2 b^3} & \frac{-a b+b^2}{a^3-3 a b^2+2 b^3} \\ \frac{-a b+b^2}{a^3-3 a b^2+2 b^3} & \frac{-a b+b^2}{a^3-3 a b^2+2 b^3} & \frac{a^2-b^2}{a^3-3 a b^2+2 b^3} \\ \end{array} \right]$$

$$\left[ \begin{array}{cccc} a & b & b & b \\ b & a & b & b \\ b & b & a & b \\ b & b & b & a \\ \end{array} \right]^{-1} = \left[ \begin{array}{cccc} \frac{a^3-3 a b^2+2 b^3}{a^4-6 a^2 b^2+8 a b^3-3 b^4} & \frac{-a^2 b+2 a b^2-b^3}{a^4-6 a^2 b^2+8 a b^3-3 b^4} & \frac{-a^2 b+2 a b^2-b^3}{a^4-6 a^2 b^2+8 a b^3-3 b^4} & \frac{-a^2 b+2 a b^2-b^3}{a^4-6 a^2 b^2+8 a b^3-3 b^4} \\ \frac{-a^2 b+2 a b^2-b^3}{a^4-6 a^2 b^2+8 a b^3-3 b^4} & \frac{a^3-3 a b^2+2 b^3}{a^4-6 a^2 b^2+8 a b^3-3 b^4} & \frac{-a^2 b+2 a b^2-b^3}{a^4-6 a^2 b^2+8 a b^3-3 b^4} & \frac{-a^2 b+2 a b^2-b^3}{a^4-6 a^2 b^2+8 a b^3-3 b^4} \\ \frac{-a^2 b+2 a b^2-b^3}{a^4-6 a^2 b^2+8 a b^3-3 b^4} & \frac{-a^2 b+2 a b^2-b^3}{a^4-6 a^2 b^2+8 a b^3-3 b^4} & \frac{a^3-3 a b^2+2 b^3}{a^4-6 a^2 b^2+8 a b^3-3 b^4} & \frac{-a^2 b+2 a b^2-b^3}{a^4-6 a^2 b^2+8 a b^3-3 b^4} \\ \frac{-a^2 b+2 a b^2-b^3}{a^4-6 a^2 b^2+8 a b^3-3 b^4} & \frac{-a^2 b+2 a b^2-b^3}{a^4-6 a^2 b^2+8 a b^3-3 b^4} & \frac{-a^2 b+2 a b^2-b^3}{a^4-6 a^2 b^2+8 a b^3-3 b^4} & \frac{a^3-3 a b^2+2 b^3}{a^4-6 a^2 b^2+8 a b^3-3 b^4} \\ \end{array} \right]$$

It appears that there should be a formula but I am not sure how to derive it. In the end, I am looking for a numerical result.

Answer 1

You can use

$$ \textbf{P} = \left[ \begin{array}{cccc} 1 & 1 & \cdots & 1\\ 1 & 1 & \cdots & 1\\ \vdots & \vdots & \ddots & \vdots\\ 1 & 1 & \cdots & 1\\ \end{array} \right] $$

Note that

$$ \textbf{P}^2 = n \textbf{P} $$

You want the inverse of

$$ b \textbf{P} + (a-b) \textbf{I} $$

You can try

$$ k \textbf{P} + \frac{1}{a-b} \textbf{I} $$

So

$$ \Big( b \textbf{P} + (a-b) \textbf{I} \Big) \Big( k \textbf{P} + \frac{1}{a-b} \textbf{I} \Big) = \textbf{I} $$

Then we get

$$ \Big( nbk + (a-b)k + \frac{b}{a-b}\Big) \textbf{P} + \textbf{I} = \textbf{I} $$

so you can solve $k$ and you find

$$ k = \frac{-b}{(a-b)(nb+a-b)} $$

So you would get

$$ \Big( b \textbf{P} + (a-b) \textbf{I} \Big)^{-1} = \frac{-b}{(a-b)(nb+a-b)} \textbf{P} + \frac{1}{a-b} \textbf{I} {}{}{} $$

Answer 2

I believe you can obtain the inverse using the Sherman-Morrison formula,

namely

$(A+uv^T)^{-1} = A^{-1} - {A^{-1}uv^T A^{-1} \over 1 + v^T A^{-1}u}$.

[See for example http://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula ]

Choosing

$$ u = v = \sqrt{b}[\begin{array} 1 & 1 & 1 & \cdots & 1 \end{array}]^T, $$ and $$ A = (a-b)I $$
yields a matrix of the correct form, i.e. $$ (A+uv^T) = \left[\begin{array}{cccc} a & b & \cdots & b\\ b & a & & \vdots\\ \vdots & &\ddots &b\\ b & \cdots & b & a \end{array}\right]. $$

The inverse is then obtained from a straightforward application of the formula.

The nice thing about obtaining the inverse using the Sherman Morrison formula is that it generalizes to an arbitrary invertible matrix $A$, not just a multiple of the identity. The only restriction is that $1+v^TA^{-1}u \neq 0$.

Answer 3

Of course, the result is well-known. Yet, using a software as mathematica, we can easily guess the correct formula. To do that, you must ask mathematica to factorize. For instance, let $n=4$ ; if $i\not=j$, then $factor((A^{-1})_{i,j})=\dfrac{-b}{(a-b)(3b+a)}$ and $factor((A^{-1})_{i,i})=\dfrac{2b+a}{(a-b)(3b+a)}$. With few values of $n$, I obtain (unless I am not good for mathematics) ${A^{-1}}_{i,j}=\dfrac{-b}{(a-b)(nb+a-b)}$ and ${A^{-1}}_{i,i}=\dfrac{nb+a-2b}{(a-b)(nb+a-b)}$. Then, it is necessary to prove that this result is correct ; yet, the work is easier when we know (before) the formula.