# ILP program to find a centrosymmetric Hadamard matrix

A question in mathoverflow asks if there exists a centrosymmetric Hadamard matrix of order 36.

An $$n \times n$$ matrix $$A = (a_{i,j})$$ is centrosymmetric if: $$a_{i,j} = a_{n-i+1, n-j+1}, \space i=1,\dots,n, j=1,\dots,n \tag0\label0$$

In a Hadamard matrix, each pair of rows has matching entries in exactly half of their columns and mismatched entries in the remaining columns, and each pair of columns has matching entries in exactly half of their rows and mismatched entries in the remaining rows. It is sufficient to enforce the property for rows or columns to have it verified for columns or rows respectively.

Using some answers from this site, I have modeled the problem as an ILP with all binary variables:

$$\min a_{1,1} \tag1\label1$$

subject to:

$$x_{k,i,j} \le a_{k,i} + a_{k,j} \le 2-x_{k,i,j}, \space 1 \le k \le n, 1 \le i \lt j \le n \tag2\label2$$ $$a_{k,i} - a_{k,j} - x_{k,i,j} \le 0, \space 1 \le k \le n, 1 \le i \lt j \le n \tag3\label3$$ $$a_{k,j} - a_{k,i} - x_{k,i,j} \le 0, \space 1 \le k \le n, 1 \le i \lt j \le n \tag4\label4$$ $$\sum_{k=1}^n x_{k,i,j} = n/2, \space 1 \le k \le n, 1 \le i \lt j \le n \tag5\label5$$

where $$\eqref{2}$$, $$\eqref{3}$$, $$\eqref{4}$$ are equivalent to $$x_{k,i,j} = 1 \iff a_{k,i} \not = a_{k,j}$$, and $$\eqref{5}$$ requires that exactly half of the entries of the two columns $$i$$ and $$j$$ match at the same row.

Values $$-1$$ in the Hadamard matrix are mapped to $$0$$ in the program formulation.

In the actual implementation, I have replaced the variables when equal due to $$\eqref{0}$$, therefore the actual variables in the implementation are the entries of the first full $$n/2$$ rows of the matrix.

An example in LP format for $$n=4$$ is:

Minimize
obj: a1,1
Subject To
c1: a1,1 + a1,2 - x1,1,2 >= 0
c2: a1,1 + a1,2 + x1,1,2 <= 2
c3: a1,1 - a1,2 - x1,1,2 <= 0
c4: a1,2 - a1,1 - x1,1,2 <= 0
c5: a2,1 + a2,2 - x2,1,2 >= 0
c6: a2,1 + a2,2 + x2,1,2 <= 2
c7: a2,1 - a2,2 - x2,1,2 <= 0
c8: a2,2 - a2,1 - x2,1,2 <= 0
c9: a2,4 + a2,3 - x3,1,2 >= 0
c10: a2,4 + a2,3 + x3,1,2 <= 2
c11: a2,4 - a2,3 - x3,1,2 <= 0
c12: a2,3 - a2,4 - x3,1,2 <= 0
c13: a1,4 + a1,3 - x4,1,2 >= 0
c14: a1,4 + a1,3 + x4,1,2 <= 2
c15: a1,4 - a1,3 - x4,1,2 <= 0
c16: a1,3 - a1,4 - x4,1,2 <= 0
c17: a1,1 + a1,3 - x1,1,3 >= 0
c18: a1,1 + a1,3 + x1,1,3 <= 2
c19: a1,1 - a1,3 - x1,1,3 <= 0
c20: a1,3 - a1,1 - x1,1,3 <= 0
c21: a2,1 + a2,3 - x2,1,3 >= 0
c22: a2,1 + a2,3 + x2,1,3 <= 2
c23: a2,1 - a2,3 - x2,1,3 <= 0
c24: a2,3 - a2,1 - x2,1,3 <= 0
c25: a2,4 + a2,2 - x3,1,3 >= 0
c26: a2,4 + a2,2 + x3,1,3 <= 2
c27: a2,4 - a2,2 - x3,1,3 <= 0
c28: a2,2 - a2,4 - x3,1,3 <= 0
c29: a1,4 + a1,2 - x4,1,3 >= 0
c30: a1,4 + a1,2 + x4,1,3 <= 2
c31: a1,4 - a1,2 - x4,1,3 <= 0
c32: a1,2 - a1,4 - x4,1,3 <= 0
c33: a1,1 + a1,4 - x1,1,4 >= 0
c34: a1,1 + a1,4 + x1,1,4 <= 2
c35: a1,1 - a1,4 - x1,1,4 <= 0
c36: a1,4 - a1,1 - x1,1,4 <= 0
c37: a2,1 + a2,4 - x2,1,4 >= 0
c38: a2,1 + a2,4 + x2,1,4 <= 2
c39: a2,1 - a2,4 - x2,1,4 <= 0
c40: a2,4 - a2,1 - x2,1,4 <= 0
c41: a2,4 + a2,1 - x3,1,4 >= 0
c42: a2,4 + a2,1 + x3,1,4 <= 2
c43: a2,4 - a2,1 - x3,1,4 <= 0
c44: a2,1 - a2,4 - x3,1,4 <= 0
c45: a1,4 + a1,1 - x4,1,4 >= 0
c46: a1,4 + a1,1 + x4,1,4 <= 2
c47: a1,4 - a1,1 - x4,1,4 <= 0
c48: a1,1 - a1,4 - x4,1,4 <= 0
c49: a1,2 + a1,3 - x1,2,3 >= 0
c50: a1,2 + a1,3 + x1,2,3 <= 2
c51: a1,2 - a1,3 - x1,2,3 <= 0
c52: a1,3 - a1,2 - x1,2,3 <= 0
c53: a2,2 + a2,3 - x2,2,3 >= 0
c54: a2,2 + a2,3 + x2,2,3 <= 2
c55: a2,2 - a2,3 - x2,2,3 <= 0
c56: a2,3 - a2,2 - x2,2,3 <= 0
c57: a2,3 + a2,2 - x3,2,3 >= 0
c58: a2,3 + a2,2 + x3,2,3 <= 2
c59: a2,3 - a2,2 - x3,2,3 <= 0
c60: a2,2 - a2,3 - x3,2,3 <= 0
c61: a1,3 + a1,2 - x4,2,3 >= 0
c62: a1,3 + a1,2 + x4,2,3 <= 2
c63: a1,3 - a1,2 - x4,2,3 <= 0
c64: a1,2 - a1,3 - x4,2,3 <= 0
c65: a1,2 + a1,4 - x1,2,4 >= 0
c66: a1,2 + a1,4 + x1,2,4 <= 2
c67: a1,2 - a1,4 - x1,2,4 <= 0
c68: a1,4 - a1,2 - x1,2,4 <= 0
c69: a2,2 + a2,4 - x2,2,4 >= 0
c70: a2,2 + a2,4 + x2,2,4 <= 2
c71: a2,2 - a2,4 - x2,2,4 <= 0
c72: a2,4 - a2,2 - x2,2,4 <= 0
c73: a2,3 + a2,1 - x3,2,4 >= 0
c74: a2,3 + a2,1 + x3,2,4 <= 2
c75: a2,3 - a2,1 - x3,2,4 <= 0
c76: a2,1 - a2,3 - x3,2,4 <= 0
c77: a1,3 + a1,1 - x4,2,4 >= 0
c78: a1,3 + a1,1 + x4,2,4 <= 2
c79: a1,3 - a1,1 - x4,2,4 <= 0
c80: a1,1 - a1,3 - x4,2,4 <= 0
c81: a1,3 + a1,4 - x1,3,4 >= 0
c82: a1,3 + a1,4 + x1,3,4 <= 2
c83: a1,3 - a1,4 - x1,3,4 <= 0
c84: a1,4 - a1,3 - x1,3,4 <= 0
c85: a2,3 + a2,4 - x2,3,4 >= 0
c86: a2,3 + a2,4 + x2,3,4 <= 2
c87: a2,3 - a2,4 - x2,3,4 <= 0
c88: a2,4 - a2,3 - x2,3,4 <= 0
c89: a2,2 + a2,1 - x3,3,4 >= 0
c90: a2,2 + a2,1 + x3,3,4 <= 2
c91: a2,2 - a2,1 - x3,3,4 <= 0
c92: a2,1 - a2,2 - x3,3,4 <= 0
c93: a1,2 + a1,1 - x4,3,4 >= 0
c94: a1,2 + a1,1 + x4,3,4 <= 2
c95: a1,2 - a1,1 - x4,3,4 <= 0
c96: a1,1 - a1,2 - x4,3,4 <= 0
c97: x1,1,2 + x2,1,2 + x3,1,2 + x4,1,2 = 2
c98: x1,1,3 + x2,1,3 + x3,1,3 + x4,1,3 = 2
c99: x1,1,4 + x2,1,4 + x3,1,4 + x4,1,4 = 2
c100: x1,2,3 + x2,2,3 + x3,2,3 + x4,2,3 = 2
c101: x1,2,4 + x2,2,4 + x3,2,4 + x4,2,4 = 2
c102: x1,3,4 + x2,3,4 + x3,3,4 + x4,3,4 = 2
Binary
a1,1 a1,2 a1,3 a1,4 a2,1 a2,2 a2,3 a2,4 x1,1,2 x2,1,2 x3,1,2 x4,1,2 x1,1,3 x2,1,3 x3,1,3 x4,1,3 x1,1,4 x2,1,4 x3,1,4 x4,1,4 x1,2,3 x2,2,3 x3,2,3 x4,2,3 x1,2,4 x2,2,4 x3,2,4 x4,2,4 x1,3,4 x2,3,4 x3,3,4 x4,3,4
End

I have uploaded here the C code for generating the ILP program. It can be run there setting as inline parameter the value of $$n$$.

I have tried some $$n=4k$$ cases successfully with up to few seconds of solving time for $$n=4,8,16$$ with Gurobi and CPLEX at NEOS Server. Case $$n=20$$ has been solved in 10 seconds by Gurobi. Strangely, CPLEX took about two hours for $$n=20$$, while it solved $$n=24$$ in a few minutes. Case $$n=24$$ with Gurobi was solved in about five hours. Cases $$n=32$$ and $$n=36$$ (Gurobi and CPLEX) timed out at NEOS after the 8 hours maximum.

Note that due to the linked question, $$n = 2^k, k \gt 1$$, is guaranteed to have a solution. The same question conjectures that $$n=4(4k+3), k \ge 0$$ has no solution.

Is there any way to optimize or parallelize the above problem to get a solution for $$n=36$$, assuming it exists? Or at least for $$n=32$$, where it surely exists?

• Rather than linearizing, you might instead try imposing quadratic constraints $\sum_k a_{ki} a_{kj} = n/2$ in place of $(2)$ through $(5)$. Nov 29, 2022 at 22:22
• Changing from $\{-1,1\}$ to $\{0,1\}$ does not preserve this property: "It is sufficient to enforce the property for rows or columns to have it verified for columns or rows respectively." For example, consider $$\begin{pmatrix}0 &0 & 0 &0 \\ 1 &1 & 1 & 1 \\ 1 &1 & 1 & 1 \\ 0 &0 & 0 &0\end{pmatrix}$$ Two workarounds: 1. Enforce the constraints for both row pairs and column pairs, 2. Replace the quadratic constraints with $\sum_k (2a_{ki}-1)(2a_{kj}-1) = 0$. Nov 30, 2022 at 1:33
• Also, you might have better luck using a constant $0$ objective so the first feasible solution found is optimal. Nov 30, 2022 at 1:38
• @RobPratt in your counterexample neither rows nor columns satisfy the requirement, because we have all elements equal or all different. I don't think it changes anything if we replace $-1$ with $0$ and just compare elements (multiplication is another thing, of course). I have taken for granted the equivalence of checking rows or columns from the Wikipedia and other sources, and verified it works in the found solutions for $n=4$ and $n=8$. Nov 30, 2022 at 7:25
• Integer programming may not be the best way to tackle the problem. While you can (and did) formulate it as an IP, the structure does not make use of the "bound" part of "branch and bound", and I'm not sure how useful the cuts produced by IP solvers would be. You might have better luck with a constraint programming model ... or it may be your best bet is just to code an enumerative search algorithm from scratch.
– prubin
Nov 30, 2022 at 22:56

In case you don't want to linearize as @robpratt suggests another way is: Initialize $$x_{kij}=0$$.
C1 = $$a_{kj}-a_{ki} \le x_{kij}$$
C2 = $$a_{ki} - a_{kj} \le x_{kij}$$.
Only thing remains when $$a_{ki} - a_{kj}=0$$, then $$x$$ is free. But if initialized to 0 and last constraint on how many $$x$$ to turn 1, solver will not change any $$x$$ to 1 unless forced.
• Because $(5)$ is an equality, the linearization must also enforce $(2)$, not just $(3)$ and $(4)$. Nov 30, 2022 at 0:23