I have a special $NM \times NM$ matrix of the form
\begin{align*}
S = \left[
\begin{array}{cccc}
V + \lambda I & V & \cdots & V \\
V & V + \lambda I & \cdots & V \\
\cdots & \cdots & \cdots & \cdots \\
V & V & \cdots & V + \lambda I
\end{array}
\right]
\end{align*}
where $V$ is a symmetric $N\times N$ matrix and $I$ is an identity matrix of size $N$.
I want to know if there is any way to express $S^{-1}$ in a simpler form involving $V^{-1}$ and $(V+\lambda I)^{-1}$.
Best Answer
Using the results from your other question and the idea of Jean Marie's answer: $S = (I_M\otimes\lambda I_N) + ({\bf 1}{\bf 1}^T \otimes V)$
We can rewrite @Accumulation's answer in the stringent form:
\begin{align} S^{-1} &= \tfrac{1}{\lambda}\big(I_{mn} - (-\tfrac{1}{\lambda}{\bf 1}_m{\bf 1}_m^T \otimes V)\big)^{-1} \\&=\tfrac{1}{\lambda}\sum_{k=0}^{\infty} (-\tfrac{1}{\lambda}{\bf 1}_m{\bf 1}_m^T \otimes V)^k \\ &=\tfrac{1}{\lambda}(I_m\otimes I_n) + \tfrac{1}{m\lambda}\sum_{k=1}^{\infty}(-\tfrac{m}{\lambda})^k ({\bf 1}_m{\bf 1}_m^T \otimes V^k) \\ &=\tfrac{1}{\lambda}(I_m\otimes I_n) + \tfrac{1}{\lambda}\big(\tfrac{1}{m}{\bf 1}_m{\bf 1}_m^T \otimes \sum_{k=1}^{\infty}(-\tfrac{m}{\lambda}V)^k\big) \\ &=\tfrac{1}{\lambda}(I_m\otimes I_n) +\tfrac{1}{\lambda}\big(\tfrac{1}{m}{\bf 1}_m{\bf 1}_m^T \otimes (I_n+\tfrac{m}{\lambda}V)^{-1}-I_n\big)\\ &= \tfrac{1}{\lambda}\Big((C_m\otimes I_n)+ \big(\tfrac{1}{m}{\bf 1}_m{\bf 1}_m^T \otimes (I_n+\tfrac{m}{\lambda}V)^{-1}\big)\Big)\end{align}
Note that the matrix $C_m = I_m - \tfrac{1}{m}{\bf 1}_m {\bf 1}_m^T $ is known as the (geometric) centering matrix
Here the Neumann series converge if $\| \tfrac{1}{\lambda}{\bf 1}_M{\bf 1}_M^T \otimes V\|<1 $ and $\|\tfrac{M}{\lambda} V \|<1$. If we take the Frobenius norm, both are equivalent to $\|V\|_{F}<\tfrac{|\lambda|}{M}$
EDIT: fixed mistake with handling $k=0$ term