# Mixed integer programming formulation for Ising model [closed]

I want to implement a minimisation on a 2D spin Ising model with 30x30 grid. The spin variables is 0,1 and the objective is to minimize the sum of products of spins. For simplicity, I only include NN pair interaction and a "triangular" interaction term. The formulation I give is as followed. And I wish to ask for suggestions whether there would be some stronger formulation that could make computation faster. Thank you :D

the original objective is:

$$\sum\limits_{i = 1}^{29} {\sum\limits_{j = 1}^{30} {{J_1}{s_{i,j}}{s_{i + 1,j}}} } + \sum\limits_{i = 1}^{30} {\sum\limits_{j = 1}^{29} {{J_1}{s_{i,j}}{s_{i,j + 1}}} } + \sum\limits_{i = 1}^{29} {\sum\limits_{j = 1}^{29} {{J_3}{s_{i,j}}{s_{i + 1,j}}{s_{i,j + 1}}} }$$

$${s_{i,j}} \in \{ 0,1\}$$

we linearize it so that:

[\begin{align} & \min : \\ & \sum\limits_{i=1}^{29}{\sum\limits_{j=1}^{30}{{{J}_{1}}{{s}_{i,j,1}}}}+\sum\limits_{i=1}^{30}{\sum\limits_{j=1}^{29}{{{J}_{1}}{{s}_{i,j,2}}}}+\sum\limits_{i=1}^{29}{\sum\limits_{j=1}^{29}{{{J}_{3}}{{s}_{i,j,3}}}} \\ & \text{subject to:} \\ & \forall i=1...29\,j=1...30\,\ k=1...4 \\ & \left[ \begin{matrix} {{s}_{i,j}} \\ {{s}_{i+1,j}} \\ {{s}_{i,j,1}} \\ \end{matrix} \right]={{y}_{i,j,1,1}}\left[ \begin{matrix} 0 \\ 0 \\ 0 \\ \end{matrix} \right]+{{y}_{i,j,1,2}}\left[ \begin{matrix} 1 \\ 0 \\ 0 \\ \end{matrix} \right]+{{y}_{i,j,1,3}}\left[ \begin{matrix} 0 \\ 1 \\ 0 \\ \end{matrix} \right]+{{y}_{i,j,1,4}}\left[ \begin{matrix} 1 \\ 1 \\ 1 \\ \end{matrix} \right] \\ & \sum\limits_{k=1}^{4}{{{y}_{i,j,1,k}}=1} \\ & 0\le y\le 1 \\ & {{y}_{i,j,1,k}}\text{ is integral} \\ & \\ & \forall i=1...30\,j=1...29\,\,k=1...4 \\ & \left[ \begin{matrix} {{s}_{i,j}} \\ {{s}_{i,j+1}} \\ {{s}_{i,j,1}} \\ \end{matrix} \right]={{y}_{i,j,2,1}}\left[ \begin{matrix} 0 \\ 0 \\ 0 \\ \end{matrix} \right]+{{y}_{i,j,2,2}}\left[ \begin{matrix} 1 \\ 0 \\ 0 \\ \end{matrix} \right]+{{y}_{i,j,2,3}}\left[ \begin{matrix} 0 \\ 1 \\ 0 \\ \end{matrix} \right]+{{y}_{i,j,2,4}}\left[ \begin{matrix} 1 \\ 1 \\ 1 \\ \end{matrix} \right] \\ & \sum\limits_{k=1}^{4}{{{y}_{i,j,2,k}}=1} \\ & 0\le {{y}_{i,j,2,k}}\le 1 \\ & {{y}_{i,j,2,k}}\text{ is integral} \\ & \\ & \forall i=1...29\,j=1...29\,\,k=1...8 \\ & \left[ \begin{matrix} {{s}_{i,j}} \\ {{s}_{i+1,j}} \\ {{s}_{i,j+1}} \\ {{s}_{i,j,3}} \\ \end{matrix} \right]={{y}_{i,j,3,1}}\left[ \begin{matrix} 0 \\ 0 \\ 0 \\ 0 \\ \end{matrix} \right]+{{y}_{i,j,3,2}}\left[ \begin{matrix} 1 \\ 0 \\ 0 \\ 0 \\ \end{matrix} \right]+{{y}_{i,j,3,3}}\left[ \begin{matrix} 0 \\ 1 \\ 0 \\ 0 \\ \end{matrix} \right]+{{y}_{i,j,3,4}}\left[ \begin{matrix} 1 \\ 1 \\ 0 \\ 0 \\ \end{matrix} \right]+{{y}_{i,j,3,5}}\left[ \begin{matrix} 0 \\ 0 \\ 1 \\ 0 \\ \end{matrix} \right]+{{y}_{i,j,3,6}}\left[ \begin{matrix} 1 \\ 0 \\ 1 \\ 0 \\ \end{matrix} \right]+{{y}_{i,j,3,7}}\left[ \begin{matrix} 0 \\ 1 \\ 1 \\ 0 \\ \end{matrix} \right]+{{y}_{i,j,3,8}}\left[ \begin{matrix} 1 \\ 1 \\ 1 \\ 1 \\ \end{matrix} \right] \\ & \sum\limits_{k=1}^{8}{{{y}_{i,j,3,k}}=1} \\ & 0\le {{y}_{i,j,3,k}}\le 1 \\ & {{y}_{i,j,3,k}}\text{ is integral} \\ & \\ & \forall i=1...30\,j=1...30 \\ & 0\le {{s}_{i,j}}\le 1 \\ & {{s}_{i,j}}\quad \text{integral} \\ \end{align}]

• You probably don't want to reinvent the wheel, especially when it comes to Ising model simulation. I would suggest looking at Metropolis (Monte Carlo) implementations or cluster algorithms like Swendsen/Wang. – Alex R. Jul 16 '14 at 16:48
• Actually, we don't want to do monte carlo on spin systems. The problem for monte carlo is that, it could not give me the true ground state for certain and it could not provide an bound on how low energy the true ground state would possibly be. Thank you :D – user40780 Jul 16 '14 at 20:51
• scicomp.stackexchange.com/q/14130/4274. Please do not post the same question on multiple sites. Each community should have an honest shot at answering without anybody's time being wasted. – D.W. Apr 28 '16 at 5:10
• I'm voting to close this question as off-topic because it was cross-posted and answered at another site. – Todd Trimble Apr 28 '16 at 11:38

• They optimize over wavefunction, which is a linear combination of all possible states, $2^{30 \times 30}$ in your case. You seem to restrict this to pure states only. – Dmitry Savostyanov Jul 21 '14 at 23:55
Actually, one way I realize is to use this kind of constraints: $${s_i} + {s_j} - 1 \le {s_i}*{s_j} \le {s_i}$$