Inverse of the matrix

Note that $$U(a,b)U(c,d) = U(ac+3bd, bc+ad+2bd)$$ This can be seen because, if we let $A$ be the matrix each of whose entries are $1$'s, since $U(a,b)=(a-b)I+bA$, and since $A^2=4A$ $$((a-b)I+bA)((c-d)I+dA)$$ $$ = (a-b)(c-d)I+(b(c-d)+d(a-b)+4bd)A$$ $$= ((ac+3bd)-(bc+ad+2bd))I+(bc+ad+2bd)A$$ So, if $U(a,b)^{-1}$ exists and is equal to $U(c,d)$ for some $c,d$, then we would have $ac+3bd=1$ and $bc+ad+2bd=0$. The second equation gives us $$bc = -(a+2b)d$$ The determinate of $U(a,b)$ is actually $(a-b)^3(a+3b)$, so suppose $U(a,b)\neq 0$, so that neither of these factors are zero. Plugging the above into the first equation, $$-a(a+2b)d+3b^2d=b$$ $$\implies d = \frac{b}{-a^2-2ab+3b^2} = -\frac{b}{(a-b)(a+3b)}$$ Finally, $$c = \frac{a+2b}{(a-b)(a+3b)}$$ Therefore, assuming $a\neq b$ and $a\neq -3b$, $$U(a,b)^{-1}=U\left(\frac{a+2b}{(a-b)(a+3b)} , -\frac{b}{(a-b)(a+3b)}\right)$$

this is $$ (a-b)I + b J, $$ where $I$ is the identity matrix and $J$ is the matrix of all ones.

The main observation is that $J^2 = n J,$ with square $n$ by $n$ matrices. Here, $$ J^2 = 4 J. $$

So, solve $$ \left( (a-b)I + b J\right) \left( xI + y J\right) = I $$ Which becomes the system of equations in $x,y$ $$ (a-b)x = 1, \; \; \; (a-b) y + bx + 4by = 0 $$ $$ x = \frac{1}{a-b}, \; \; \; (a+3b) y + \frac{b}{a-b} = 0, $$ $$ x = \frac{1}{a-b}, \; \; \; y = \frac{-b}{(a-b)(a+3b)}. $$ THE REQUESTED INVERSE IS $$ \frac{1}{a-b} I - \frac{b}{(a-b)(a+3b)} J $$ In particular, the resulting diagonal entries (which are $x+y$) are $$ \frac{a +2b}{(a-b)(a+3b)} $$

Note that this is sometimes impossible; it is easy to find the eigenvalues of the original matrix. That is, we can tell easily when it has no inverse. The eigenvalues are, simply, $(a+3b,a-b,a-b,a-b).$

The columns of the matrix below are pairwise perpendicular and are eigenvectors of the given matrix. This fails to be an orthogonal matrix in that the columns are of varying lengths, but that can be fixed by dividing the columns by four separate numbers, namely $2, \sqrt 2, \sqrt 6, \sqrt {12}$ $$ \left( \begin{array}{rrrr} 1 & -1 & -1 & -1 \\ 1 & 1 & -1 & -1 \\ 1 & 0 & 2 & -1 \\ 1 & 0 & 0 & 3 \end{array} \right). $$

Note $\, U = a + b x + b x^2 + b x^3\, $ for a circulant $x$ satisfying $\,0 = x^4 - 1 = (x\!-\!1)\ (x\!+\!1)(x^2\!+\!1)$. Thus it suffices to compute $U^{-1}\!\!\pmod{x^4\!-1}$ by CRT, working modulo the above listed factors.

${\rm mod}\ \color{#0a0}{(x\!+\!1)(x^2\!+\!1)}\!:\,\ U^{-1}\equiv \color{#c00}{\dfrac{1}{a\!-\!b}}\ $ by $\ x\equiv -1 $ or $\,x^2\equiv -1\,\Rightarrow\, U\equiv a\!-\!b\,$

${\rm mod}\ x\!-\!1\!:\, \ \dfrac{1}{a\!+\!3b}\equiv U^{-1}\equiv \color{#c00}{\dfrac{1}{a\!-\!b}} + c\:\color{#0a0}{(x\!+\!1)(x^2\!+\!1)}\equiv \dfrac{1}{a\!-\!b}+4c\ $ by $\ x\equiv 1$

hence $\ 4c \equiv \dfrac{1}{a\!+\!3b}-\dfrac{1}{a\!-\!b}\ $ therefore $\ c \equiv \dfrac{-b\ \ \ }{(a\!-\!b)(a\!+\!3b)}$

Thus $\ U^{-1} = \underbrace{\dfrac{1}{a\!-\!b} + \color{#0a0}{c}}_{\large d\,} +\color{#0a0}{c\,(x\!+\!x^2\!+\!x^3)}\, =\, \bbox[6px,border:1px solid red]{U(d,\,c)}\,,\ $ $\, d = \dfrac{a\!+\!2b}{(a\!-\!b)(a\!+\!3b)}$

I will use the structure of a circulant matrix.

Let $T$ be the matrix $$ T=\left(\begin{array}{cccc}0&1&0&0\\0&0&1&0\\0&0&0&1\\1&0&0&0 \end{array}\right). $$ We have $T^4=I$, and your matrix is $$ U(a,b)=C(T):=aI+bT+bT^2+bT^3. $$ The matrix $T$ has eigenvalues $1,i,-1,-i$, so there is a complex matrix $P$ such that $P^{-1}TP=diag(1,i,-1,-i)$. Actually it is easy to see that (the discrete Fourier transform 4x4 matrix $$ P=\left(\begin{array}{cccc} 1&1&1&1\\ 1&i&-1&-i\\ 1&-1&1&-1\\ 1&-i&-1&i \end{array}\right) $$ works. The same matrix $P$ diagonalizes all the powers $T^j$ and $U(a,b)$ has eigenvalues $C(1)=a+3b$, $C(i)=a-b$, $C(-1)=a-b$ and $C(-i)=a-b$, so $$ P^{-1}U(a,b)P=diag(a+3b,a-b,a-b,a-b). $$ From this we see that $$ U(a,b)^{-1}=P\ diag(1/(a+3b),1/(a-b),1/(a-b),1/(a-b))\ P^{-1}. $$ If you find a cubic polynomial $Q(T)$ such that $Q(1)=1/(a+3b)$ and $Q(\pm i)=Q(-1)=1/(a-b)$, then you can deduce that $$ U(a,b)^{-1}=Q(T). $$