Consider $M$ an $nN\times nN$ block matrix which can be written as $n\times n$ blocks, with all the "diagonal" blocks equal $A = a I$ and all the "off-diagonal" blocks equal $B = b I$ where $I$ is the $N$-dimensional identity matrix and $a, b > 0$:
\begin{bmatrix} A & B & \cdots \\
B & A & B &\cdots \\
\vdots & B & A &\cdots \\
& \vdots & B & \ddots & \\ \end{bmatrix}
Is there a closed form solution for the Cholesky decomposition of $M$? I guess one could iteratively use the the 2 by 2 block decomposition as outlined here https://scicomp.stackexchange.com/questions/5050/cholesky-factorization-of-block-matrices but I feel there should be a more straightforward solution.
Best Answer
Your matrix can be expressed in the form $M = P \otimes I$, where $$ P = \pmatrix{a&b&\cdots\\b&a&b&\cdots\\ \vdots&\ddots&\ddots&\ddots} $$ and $\otimes$ denotes the Kronecker product.
First, let's remark on when the matrix $M$ has a Cholesky decomposition. $M$ has a Cholesky decomposition if and only if it is positive semidefinite, and it turns out (as a consequence of the properties of the Kronecker product) that $M$ is positive semidefinite if and only if $P$ is positive semidefinite.
With that said, we'll find a Cholesky decomposition of $P$ (if $P$ is positive definite) and use this to get a Cholesky decomposition of $M$. As discussed here, it is relatively easy to show that the eigenvalues of $P$ will be equal to $a + (n-1)b$ (with multiplicity $1$) and $a-b$ (with multiplicity $n-1$). If both of these are positive, then one can compute a Choelsy decomposition for $P$.
Suppose that $L$ is a lower triangular matrix such that $LL^T = P$. It follows from the properties of the Kronecker product that $L \otimes I$ is lower triangular and that $$ (L \otimes I)(L \otimes I)^T = (LL^T) \otimes (II)^T = P \otimes I. $$ That is to say, if $P = LL^T$ is the Choelsky decomposition of $P$, then $M = (L \otimes I)(L \otimes I)^T$ will be the Cholesky decomposition of $M$.
If you are simply looking for any matrix $R$ such that $M = RR^T$, then we can avoid numerically computing the Cholesky decomposition by using the eigenvalues/eigenvectors of $P$ to produce the (symmetric) positive semidefinite square root of $M$. To start, we will find an expression for the positive semidefinite square root of $P$. Let $J$ denote the $n\times n$ matrix $$ J = \frac 1n\pmatrix{1 & 1&\cdots\\1&1&\cdots\\\vdots&\vdots&\ddots}. $$ Notably, we can express $P$ as $$ P = (a - b) I + nb J. $$ Now, write $$ Q = p I + q J. $$ Using the fact that $J^2 = J$, we find that $$ Q^2 = p^2 I + (2pq + q^2) J. $$ Thus, if we select non-negative $p,q$ such that $$ \begin{cases} p^2 = a-b\\ 2pq + q^2 = nb, \end{cases} $$ then we will find that $Q$ is positive semidefinite and that $Q^2 = P$. Solving these equations yields non-negative roots $$ p = \sqrt{a-b}, \\ q^2 + 2pq - nb = 0 \implies q = \sqrt{a-b} + \sqrt{a + (n-1)b}. $$ With that, we can show that the matrix $R = Q \otimes I$ is also positive semidefinite and satisfies $R^2 = M$.