I am trying to construct a Hadamard matrix of order 28 using Williamson's construction. But I am unable able to construct the necessary symmetric and commuting matrices.
Definition: $H_n \in M_n(1,-1)$ is a Hadamard matrix if $H_nH_n^T = nI_n $.
Theorem: If there exists symmetric commuting matrices $A,B,C,D \in M_n(1,-1)$ such that $ A^2+B^2+C^2+D^2=4nI_n $ then,
$H_{4n} =\begin{bmatrix}-A & B & C &D \\ B&A&D&-C \\C&-D&A&B \\ D&C&-B&A \end{bmatrix}$ is a Hadamard matrix.
What methods are there for constructing $A,B,C $ and $D$ ?
One idea is to restrict attention to circulant matrices, since they will automatically be commuting. This is so because a circulant matrix is a polynomial in a matrix of the form $$ \begin{bmatrix} 0 & 1 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 1 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 1\\ 1 & 0 & 0 & 0 & 0 & 0 & 0\\ \end{bmatrix}. $$
Since the matrices are also required to be symmetric, their first rows must be of the form $$ \begin{bmatrix}a & b & c & d & d & c & b\end{bmatrix} $$ in the $7\times7$ case. There are only $16$ such matrices with elements $-1$ and $1$. So the search space ends up not being very large for a computer: $16^4=65536$. This method ends up producing lots of solutions.
You can actually make this much more efficient, to the point where it can be done by hand. The squares of the row sums of the four matrices must add to $28$. This follows by the following argument. Define $$ j=\begin{bmatrix}1 & 1 & 1 & 1 & 1 & 1 & 1\end{bmatrix}. $$ Then $jA=r_Aj$ where $r_A$ is the sum of the elements of a column of $A$. (The sum is the same for all columns and is also equal to the row sums because $A$ is circulant.) Multiplying $A^2+B^2+C^2+D^2=4nI$ by $j$ on the left and by $j^T$ on the right gives $$ nr_A^2+nr_B^2+nr_C^2+nr_D^2=4n^2, $$ from which it follows that $4n$ is the sum of the squares of the four row sums.
There are only three ways to express $28$ as the sum of four squares: $$ \begin{aligned} 28&=5^2+1^2+1^2+1^2\\ &=4^2+2^2+2^2+2^2\\ &=1^2+3^2+3^2+3^2. \end{aligned} $$ The middle one doesn't work, since $7\times7$ matrices have odd row sums. Since if $(A,B,C,D)$ is a solution, then all of $(\pm A,\pm B,\pm C,\pm D)$ are solutions as well, we may assume matrices with positive row sums. Furthermore, since any permutation of the four matrices in a solution gives another solution, we may assume that either $A$ has row sum $5$ and $B$, $C$, $D$ have row sums $1$, or $A$ has row sum $1$ and $B$, $C$, $D$ have row sums $3$.
There is only one circulant symmetric $7\times7$ matrix with row sum $5$, generated by first row $$ r=\begin{bmatrix}-1 & 1 & 1 & 1 & 1 & 1 & 1\end{bmatrix}. $$ A circulant symmetric $7\times7$ matrix with row sum $1$ has four $1$s and three $-1$s in each row, and, for the $$ \begin{bmatrix}a & b & c & d & d & c & b\end{bmatrix} $$ pattern to hold, $a$ must equal $-1$, which leaves only three possible first rows, $$ \begin{align} s_1&=\begin{bmatrix}-1 & -1& 1 & 1 & 1 & 1 & -1\end{bmatrix},\\ s_2&=\begin{bmatrix}-1 & 1 & -1 & 1 & 1 & -1 & 1\end{bmatrix},\\ s_3&=\begin{bmatrix}-1 & 1 & 1 & -1 & -1 & 1 & 1\end{bmatrix}.\\ \end{align} $$ Similarly there are only three possible first rows that give row sum $3$: $$ \begin{align} t_1&=\begin{bmatrix}1 & -1& 1 & 1 & 1 & 1 & -1\end{bmatrix},\\ t_2&=\begin{bmatrix}1 & 1 & -1 & 1 & 1 & -1 & 1\end{bmatrix},\\ t_3&=\begin{bmatrix}1 & 1 & 1 & -1 & -1 & 1 & 1\end{bmatrix}.\\ \end{align} $$
If $A$ is circulant, then so is $A^2$. (It too is a polynomial in the same generating matrix as $A$ is.) So we need not do any matrix computations at all. We need only compute the first rows of $A^2$, $B^2$, $C^2$, $D^2$, or equivalently, the periodic autocorrelation functions of the first rows of $A$, $B$, $C$, and $D$, and check whether any four of them sum to $$ \begin{bmatrix}28 & 0 & 0 & 0 & 0 & 0 & 0\end{bmatrix}. $$ (The periodic autocorrelation function of the vector $\begin{bmatrix}a_0 & a_1 & \ldots & a_{n-1}\end{bmatrix}$ is defined to be the vector whose $j^\text{th}$ element, $0\le j<n$, is $$ \sum_{i=0}^{n-1}a_ia_{i+j}, $$ where the index $i+j$ is to be computed modulo $n$.)
The periodic autocorrelation function of $r$ is $$ \begin{bmatrix}7 & 3 & 3 & 3 & 3 & 3 & 3\end{bmatrix}. $$ The periodic autocorrelation functions of $s_1$, $s_2$, $s_3$ are $$ \begin{align} &\begin{bmatrix}7 & 3 & -1 & -5 & -5 & -1 & 3\end{bmatrix},\\ &\begin{bmatrix}7 & -5 & 3 & -1 & -1 & 3 & -5\end{bmatrix},\\ &\begin{bmatrix}7 & -1 & -5 & 3 & 3 & -5 & -1\end{bmatrix}.\\ \end{align} $$ The periodic autocorrelation functions of $t_1$, $t_2$, $t_3$ are $$ \begin{align} &\begin{bmatrix}7 & -1 & 3 & -1 & -1 & 3 & -1\end{bmatrix},\\ &\begin{bmatrix}7 & -1 & -1 & 3 & 3 & -1 & -1\end{bmatrix},\\ &\begin{bmatrix}7 & 3 & -1 & -1 & -1 & -1 & 3\end{bmatrix}.\\ \end{align} $$
One sees that we get solutions by taking the generating rows of $A$, $B$, $C$, $D$ to be one of the following: $$ \begin{aligned} &r,\ s_1,\ s_2,\ s_3,& &\text{with row sums 5, 1, 1, 1}\\ &s_1,\ t_1,\ t_2,\ t_2,& &\text{with row sums 1, 3, 3, 3}\\ &s_2,\ t_2,\ t_3,\ t_3,& &\text{with row sums 1, 3, 3, 3}\\ &s_3,\ t_3,\ t_1,\ t_1,& &\text{with row sums 1, 3, 3, 3}. \end{aligned} $$