Linearizing multiple XOR operations

I have never encountered this in a model. but nevertheless this is an interesting question.

How can we handle \[y = x_1 \oplus x_2 \oplus \cdots \oplus x_n\] in a linear MIP model? Here \(\oplus\) indicates the xor operation. 

Well, when we look at \(x_1 \oplus x_2 \oplus \cdots \oplus x_n\) we can see that this is identical to \(\sum x_i\) being odd. A small example illustrates that:

Checking if an expression is odd can be done by seeing if it can not be expressed as \(2z\) where \(z\) is an integer variable. So combining this, we can do:

Linearization of \(y = x_1 \oplus x_2 \oplus \cdots \oplus x_n\)
\[\begin{align} & \color{DarkRed}y = \sum_i \color{DarkRed}x_i - 2 \color{DarkRed}z \\ & \color{DarkRed}y, \color{DarkRed} x_i \in\{0,1\} \\ & \color{DarkRed}z \in \{0,1,2,\dots\} \end{align}\]

If we want, we can relax \(y\) to be continuous between 0 and 1. An upper bound on z can be
\[z \le \left \lfloor{ \frac{n}{2} }\right \rfloor \].

References

1. Express XOR with multiple inputs in zero-one integer linear programming (ILP), https://cs.stackexchange.com/questions/40737/express-xor-with-multiple-inputs-in-zero-one-integer-linear-programming-ilp
2. XOR as linear inequalities, https://yetanothermathprogrammingconsultant.blogspot.com/2016/02/xor-as-linear-inequalities.html
3. A variant of the Lights Out game, https://yetanothermathprogrammingconsultant.blogspot.com/2016/01/a-variant-of-lights-out-game.html

An interesting linearization involving xor

XOR gate sold by Amazon [1]

In [2] a question was asked:

How can we solve \[x^TAy = b\] for \(x_i, y_i \in \{-1,+1\}\)?

In this post I will discuss three different formulations:

1. A simple quadratic formulation yielding a non-convex MIQCP model
2. A MIP model based on a standard linearization
3. A MIP model based on a linearization using an xor construct

Test Data

To help with some test models we need some data for \(A\) and \(b\). Here is an instance that actually has a solution (printed with 3 decimals):

----     13 PARAMETER a  

            i1          i2          i3          i4          i5

i1       0.998       0.579       0.991       0.762       0.131
i2       0.640       0.160       0.250       0.669       0.435
i3       0.360       0.351       0.131       0.150       0.589
i4       0.831       0.231       0.666       0.776       0.304
i5       0.110       0.502       0.160       0.872       0.265


----     13 PARAMETER b                    =        2.222 

Notice that solutions are not unique: if we have a solution \((x,y)\) then another solution is \((-x,-y)\).

MIQCP Model

The simplest is to use a Mixed-Integer Quadratically Constrained (MIQCP) model that handles the quadratic equation directly. Of course, in optimization we prefer binary variables \(z \in\{0,1\}\) instead of \(z \in\{-1,+1\}\). That is easily fixed: introduce binary variables \(p_i, q_i\) and write: \[\begin{align}&x_i = 2p_i-1\\&y_i = 2q_i-1\\&p_i,q_i \in\{0,1\}\end{align}\] This maps a binary variable to \(\{-1,+1\}\). So a MIQCP model can look like:

MIQCP Formulation
\[\begin{align} & \sum_{i,j} \color{DarkRed} x_i \color{DarkRed}y_j \color{DarkBlue}a_{i,j} = \color{DarkBlue}b \\ & \color{DarkRed}x_{i} = 2\color{DarkRed}p_{i}-1\\ & \color{DarkRed}y_{i} = 2\color{DarkRed}q_{i}-1\\ & \color{DarkRed}p_i, \color{DarkRed}q_i \in\{0,1\} \end{align}\]

This is a non-convex problem, so you need a global solver like Baron, Couenne or Antigone. Most quadratic solvers only support the much easier, convex case. When expressing this feasibility problem as an optimization problem, we need to add a dummy objective. The solution can look like:

----     55 VARIABLE p.L  

i1 1.000,    i4 1.000,    i5 1.000


----     55 VARIABLE q.L  

i1 1.000,    i2 1.000,    i4 1.000


----     55 VARIABLE x.L  

i1  1.000,    i2 -1.000,    i3 -1.000,    i4  1.000,    i5  1.000


----     55 VARIABLE y.L  

i1  1.000,    i2  1.000,    i3 -1.000,    i4  1.000,    i5 -1.000


----     59 PARAMETER b2                   =        2.222  x'Ay using solution values

We find a solution that when plugged into \(x^TAy\), reproduces our right-hand side.

Linearization 1

This problem can be linearized. This will allow us to use a linear MIP solver. One way is to write:\[\sum_{i,j} x_i y_j a_{i,j} = \sum_{i,j} (2p_i-1)(2q_j-1) a_{i,j} = \sum_{i,j} \left( 4 p_i q_j - 2 p_i -2q_j + 1 \right) a_{i,j} \]

The binary product \(v_{i,j} = p_i q_j\) can be linearized using a standard formulation [3]: \[\begin{align} &v_{i,j} \le p_i\\ &v_{i,j} \le q_j \\ & v_{i,j} \ge p_i + q_j -1 \\ & p_i, q_j, v_{i,j} \in \{0,1\} \end{align}\] Combining this leads to:

MIP Formulation I

\[\begin{align} & \sum_{i,j} \left( 4 \color{DarkRed} v_{i,j} - 2 \color{DarkRed}p_i -2 \color{DarkRed} q_j + 1 \right) \color{DarkBlue} a_{i,j} = \color{DarkBlue}b \\ & \color{DarkRed}v_{i,j} \le \color{DarkRed}p_i\\ & \color{DarkRed}v_{i,j} \le \color{DarkRed}q_j \\ & \color{DarkRed}v_{i,j} \ge \color{DarkRed}p_i + \color{DarkRed}q_j -1 \\ & \color{DarkRed}x_{i} = 2\color{DarkRed}p_{i}-1\\ & \color{DarkRed}y_{i} = 2\color{DarkRed}q_{i}-1\\ & \color{DarkRed}p_i, \color{DarkRed}q_i, \color{DarkRed}v_{i,j} \in\{0,1\} \end{align}\]

It is noted that \(v\) can be relaxed to be continuous between 0 and 1: \(v_{i,j} \in [0,1]\). The solution can look like:

----     62 VARIABLE p.L  

i1 1.000,    i4 1.000,    i5 1.000


----     62 VARIABLE q.L  

i1 1.000,    i2 1.000,    i4 1.000


----     62 VARIABLE v.L  v(i,j)=p(i)*q(j)

            i1          i2          i4

i1       1.000       1.000       1.000
i4       1.000       1.000       1.000
i5       1.000       1.000       1.000


----     62 VARIABLE x.L  

i1  1.000,    i2 -1.000,    i3 -1.000,    i4  1.000,    i5  1.000


----     62 VARIABLE y.L  

i1  1.000,    i2  1.000,    i3 -1.000,    i4  1.000,    i5 -1.000


----     66 PARAMETER b2                   =        2.222  x'Ay using solution values

Note that zeroes are not printed here (so the matrix \(v\) may look a bit strange).

Linearization 2

There is another linearization which is more interesting. For this we will use the binary xor operator.

There are different ways xor is denoted. I believe the most common ones are \[\begin{align} & z = x \textbf{ xor } y\\ & z = \textbf{xor}(x,y) \\ & z = x \oplus y \\ & z = x \veebar y\end{align}\] The xor truth table is: \[\begin{matrix} \textbf{ x } & \textbf{ y } & \textbf{ z }\\ 0 & 0 & 0 \\ 0 & 1 & 1 \\ 1 & 0 & 1\\ 1 & 1 & 0
\end{matrix}\]I.e. \(z = x \textbf{ xor } y \) is 1 if and only if \(x\) and \(y\) are different.

Now consider the following thruth-table:

----     15 PARAMETER t  truth-table

             x           y         x*y           p           q     p xor q   1 - 2 xor

i1          -1          -1           1                                               1
i2          -1           1          -1                       1           1          -1
i3           1          -1          -1           1                       1          -1
i4           1           1           1           1           1                       1

From this we can see that \[x_i  y_j = 1 - 2 (p_i \textbf{ xor } q_j) \] So our quadratic constraint becomes \[\sum_{i,j} x_i y_j a_{i,j} = \sum_{i,j} \left( 1 - 2 (p_i \textbf{ xor } q_j) \right)  a_{i,j} = b\] The operation \(w_{i,j} = p_i \textbf{ xor } q_j\)  can be linearized [4]: \[\begin{align} & w_{i,j} \le p_i + q_j \\& w_{i,j} \ge p_i - q_j \\& w_{i,j} \ge q_j - p_i \\& w_{i,j} \le 2 - p_i - q_j \\ & p_i, q_j, w_{i,j} \in \{0,1\} \end{align}\] With these tools we can formulate a different linearization:

MIP Formulation II

\[\begin{align} & \sum_{i,j} \left( 1 - 2 \color{DarkRed} w_{i,j} \right) \color{DarkBlue} a_{i,j} = \color{DarkBlue}b \\ & \color{DarkRed}w_{i,j} \le \color{DarkRed}p_i + \color{DarkRed}q_j\\ & \color{DarkRed}w_{i,j} \ge \color{DarkRed}p_i - \color{DarkRed}q_j \\ & \color{DarkRed}w_{i,j} \ge \color{DarkRed}q_j - \color{DarkRed}p_i \\ & \color{DarkRed}w_{i,j} \le 2 - \color{DarkRed}p_i - \color{DarkRed}q_j \\ & \color{DarkRed}x_{i} = 2\color{DarkRed}p_{i}-1\\ & \color{DarkRed}y_{i} = 2\color{DarkRed}q_{i}-1\\ & \color{DarkRed}p_i, \color{DarkRed}q_i, \color{DarkRed}w_{i,j} \in\{0,1\} \end{align}\]

It is noted that \(w\) can be relaxed to be continuous between 0 and 1: \(w_{i,j} \in [0,1]\). A solution can look like:

----     62 VARIABLE p.L  

i1 1.000,    i4 1.000,    i5 1.000


----     62 VARIABLE q.L  

i1 1.000,    i2 1.000,    i4 1.000


----     62 VARIABLE w.L  w(i,j) = p(i) xor q(j)

            i1          i2          i3          i4          i5

i1                               1.000                   1.000
i2       1.000       1.000                   1.000
i3       1.000       1.000                   1.000
i4                               1.000                   1.000
i5                               1.000                   1.000


----     62 VARIABLE x.L  

i1  1.000,    i2 -1.000,    i3 -1.000,    i4  1.000,    i5  1.000


----     62 VARIABLE y.L  

i1  1.000,    i2  1.000,    i3 -1.000,    i4  1.000,    i5 -1.000


----     66 PARAMETER b2                   =        2.222  x'Ay using solution values

Conclusion

I am somewhat surprised how different the two linearizations look. You would not easily recognize that the two different linear formulations are essentially the same. I don't remember many times I could use an xor operation in an optimization (I used it to solve Takuzu puzzles in [5]). The second linearization earns extra points for using xor!

This question was more interesting than I anticipated.

References

1. A Fairchild/ON Semiconductor XOR gate sold by Amazon for $1.69. https://www.amazon.com/Logic-Gates-2-Input-XOR-Gate/dp/B00MEKUZV8
2. Matrix Equation & Integer Programming, https://math.stackexchange.com/questions/3047086/matrix-equation-integer-programming/3047835
3. Multiplication of Binary Variables, https://yetanothermathprogrammingconsultant.blogspot.com/2008/05/multiplication-of-binary-variables.html
4. XOR as linear inequalities, https://yetanothermathprogrammingconsultant.blogspot.com/2016/02/xor-as-linear-inequalities.html
5. Solving Takuzu puzzles as a MIP using xor, https://yetanothermathprogrammingconsultant.blogspot.com/2017/01/solving-takuzu-puzzles-as-mip-using-xor.html

More queens problems (2)

In [1], I discussed the Minimum Dominating Set of Queens Problem. Here is a variant:

Given an \(n \times n\) chess board, place \(n\) queens on the board such that the number of squares not under attack is maximized. 

For an \(n=8\) problem, the maximum number of non-dominated squares is 11. An example is:

11 white pawns are not under attack

We base the model on the one in [1]. First we define two sets of binary variables:

\[\begin{align} & x_{i,j} = \begin{cases} 1 & \text{if square \((i,j)\) is occupied by a queen}\\ 0 & \text{otherwise}\end{cases} \\ & y_{i,j} = \begin{cases} 1 & \text{if square \((i,j)\) is under attack} \\ 0 & \text{otherwise}\end{cases}\end{align}\]

The model can look like:

Maximum Non-Dominated Squares Problem

\[\begin{align}\max& \>\color{DarkRed}z=\sum_{i,j} \left(1-\color{DarkRed} y_{i,j}\right)\\ & \sum_{j'} \color{DarkRed}x_{i,j'} + \sum_{i'} \color{DarkRed}x_{i',j} + \sum_{i',j'|i'-j'=i-j} \color{DarkRed}x_{i',j'} +\sum_{i',j'|i'+j'=i+j} \color{DarkRed}x_{i',j'}-3 \color{DarkRed} x_{i,j} \ge \color{DarkBlue}M \cdot \color{DarkRed}y_{i,j} && \forall i,j \\  & \sum_{i,j}  \color{DarkRed} x_{i,j} = n \\& \color{DarkRed}x_{i,j},\color{DarkRed}y_{i,j}\in\{0,1\} \end{align}\]

The funny term \(-3x_{i,j}\) is to compensate for double counting occurrences of \(x_{i,j}\) in the previous terms [1]. We need to find a reasonable value for \(M\). An upper bound is \(M=n\). For \(n=8\), there are 48 different optimal solutions.

The problem \(n=9\) has one structural optimal solution (with 18 unattacked positions), and 3 variants by symmetry:

Is there a good way to find only structural solutions? One way is to use your own cuts as in [4], If we want to use the solution pool technology, we need to invent a constraint that forbids some symmetry (or have an objective that favors some of symmetric versions).

References

1. More queens problems, https://yetanothermathprogrammingconsultant.blogspot.com/2018/12/more-queens-problems.html
2. Bernard Lemaire, Pavel Vitushinkiy, Placing \(n\) non-dominating queens on the \(n\times n\) Chessboard, 2011, https://www.ffjm.org/upload/fichiers/N_NON_DOMINATING_QUEENS.pdf
3. Maximum Queens Chess Problem, https://www.gams.com/latest/gamslib_ml/queens.103

More queens problems

In [1] the famous \(n\)-queens problem was discussed: how many queens can we place on a chess-board such that non of them are attacked by others. Interestingly, we had some non-obvious problem with Cplex enumerating the different configurations (due to tolerance issues).

A different problem is the Minimum Dominating Set of Queens Problem:

Find the minimum number of queens needed such that they dominate each square.

Here dominate means: either a square is under attack by at least one queen, or it is occupied by a queen.

To develop a Mixed Integer Programming model, we introduce the familiar binary variables:

\[x_{i,j} = \begin{cases} 1 & \text{if square \((i,j)\) is occupied by a queen}\\ 0 & \text{otherwise}\end{cases}\]

A square \((i,j)\) is dominated if:

1. \(x_{i,j}=1\): a queen is located here. This case will be covered by the following other cases, so we don't have to worry about this.
2. There is a queen in row \(i\): \[\sum_{j'} x_{i,j'} \ge 1\]
3. There is a queen in column \(j\):\[\sum_{i'} x_{i',j} \ge 1\]
4. There is a queen in the same diagonal. The diagonal is described by \((i',j')\) such that \(i'-j'=i-j\). So we have: \[\sum_{i',j'|i'-j'=i-j} x_{i',j'}\ge 1\] If you prefer you can write this as \[\sum_{i'|1 \le i'-i+j \le n } x_{i',i'-i+j} \ge 1 \] 
5. There is a queen in the same anti-diagonal. The anti-diagonal is described by \((i',j')\) such that \(i'+j'=i+j\). So we have: \[\sum_{i',j'|i'+j'=i+j} x_{i',j'}\ge 1\] This can also be written as \[\sum_{i'| 1 \le i+j-i' \le n} x_{i',i+j-i'} \ge 1\]

We can write these conditions in one big constraint. So we can write:

Minimum Dominating Set of Queens Problem
\[\begin{align}\min& \>\color{DarkRed}z=\sum_{i,j} \color{DarkRed} x_{i,j}\\ & \sum_{j'} \color{DarkRed}x_{i,j'} + \sum_{i'} \color{DarkRed}x_{i',j} + \sum_{i',j'|i'-j'=i-j} \color{DarkRed}x_{i',j'} +\sum_{i',j'|i'+j'=i+j} \color{DarkRed}x_{i',j'} \ge 1 && \forall i,j \\ & \color{DarkRed}x_{i,j}\in\{0,1\} \end{align}\]

For an \(8\times 8\) board, we need 5 queens. Here are two solutions:

You can verify that every unoccupied square is under attack.

Note that this formulation actually does a bit of double counting. E.g. when we look at the constraint for cell \((2,3)\) we see:

e(2,3)..  x(1,2) + x(1,3) + x(1,4) + x(2,1) + x(2,2) + 4*x(2,3) + x(2,4)
     
      + x(2,5) + x(2,6) + x(2,7) + x(2,8) + x(3,2) + x(3,3) + x(3,4) + x(4,1)
     
      + x(4,3) + x(4,5) + x(5,3) + x(5,6) + x(6,3) + x(6,7) + x(7,3) + x(7,8)
     
      + x(8,3) =G= 1 ; (LHS = 0, INFES = 1 ****)

We see that x(2,3) has a coefficient 4. This is harmless. However it looks a bit unpolished. With a little bit more careful modeling we can get rid of this coefficient 4. A complicated way would be: \[ \sum_{j'}x_{i,j'} + \sum_{i'|i'\ne i} x_{i',j} + \sum_{i',j'|i'-j'=i-j,i'\ne i, j'\ne j} x_{i',j'} +\sum_{i',j'|i'+j'=i+j,i'\ne i, j'\ne j} x_{i',j'} \ge 1 \>\> \forall i,j\] The first summation includes \(x_{i,j}\) while the other three summations exclude explicitly an occurrence of \(x_{i,j}\). Simpler is just to subtract \(3 x_{i,j}\): \[\sum_{j'} x_{i,j'} + \sum_{i'} x_{i',j} + \sum_{i',j'|i'-j'=i-j} x_{i',j'} +\sum_{i',j'|i'+j'=i+j} x_{i',j'}-3x_{i,j} \ge 1\>\> \forall i,j\] Here we just compensate for the double counting.

Now we see:

e(2,3)..  x(1,2) + x(1,3) + x(1,4) + x(2,1) + x(2,2) + x(2,3) + x(2,4) + x(2,5)
     
      + x(2,6) + x(2,7) + x(2,8) + x(3,2) + x(3,3) + x(3,4) + x(4,1) + x(4,3)
     
      + x(4,5) + x(5,3) + x(5,6) + x(6,3) + x(6,7) + x(7,3) + x(7,8) + x(8,3)

      =G= 1 ; (LHS = 0, INFES = 1 ****)

When we ask: how many different solution are there?, we use again the Cplex solution pool. And again, we have our tolerance problem:

pool.absgap	Number of solutions
0	2
0.5	4860

As the objective is integer valued, in some sense an absolute gap tolerance of 0.5 seems the safest. With this we find the correct number of solutions [2]. See [1] for a more in-depth discussion of this tolerance problem.

References

1. Chess and solution pool, http://yetanothermathprogrammingconsultant.blogspot.com/2018/11/chess-and-solution-pool.html.
2. Weisstein, Eric W.,  Queens Problem, http://mathworld.wolfram.com/QueensProblem.html
3. Henning Fernau, minimum dominating set of queens: A trivial programming exercise?, Discrete Applied Mathematics, Volume 158, Issue 4, 28 February 2010, Pages 308-318. This is about programming instead of modeling. The issues are very different.
4. Saad Alharbi and, Ibrahim Venkat, A Genetic Algorithm Based Approach for Solving the Minimum Dominating Set of Queens Problem, Hindawi Journal of Optimization, Volume 2017. I am not sure why a meta-heuristic is a good choice for solving this problem: you will not find proven optimal solutions.