Finding magic squares with the Douglas-Rachford algorithm

TL;DR

This paper demonstrates how the Douglas-Rachford algorithm can be effectively used as a heuristic to find magic squares, despite its convergence guarantees only in convex cases.

Contribution

It introduces two novel nonconvex formulations for magic square problems and applies the Douglas-Rachford algorithm as a heuristic approach.

Findings

Douglas-Rachford successfully finds magic squares in tested formulations

The algorithm performs well despite nonconvexity of the problem

Potential for applying the method to other nonconvex puzzles

Abstract

In this expository paper, we show how to use the Douglas-Rachford algorithm as a successful heuristic for finding magic squares. The Douglas-Rachford algorithm is an iterative projection method for solving feasibility problems. Although its convergence is only guaranteed in the convex setting, the algorithm has been successfully applied to a number of similar nonconvex problems, such as solving Sudoku puzzles. We present two formulations of the nonconvex feasibility problem of finding magic squares, which are inspired by those of Sudoku, and test the Douglas-Rachford algorithm on them.

Tables

Table 1. Table 1: Magic squares found in less than 30 30 30 minutes and mean (max) time

	Binary			Integer
Order	Success	Time		Success	Time
3	99	$0.05$	$\left(0.17\right)$	100	$0.01$	$\left(0.02\right)$
4	93	$1.67$	$\left(13.02\right)$	64	$0.06$	$\left(2.64\right)$
5	79	$436.86$	$\left(1739.93\right)$	59	$0.67$	$\left(5.63\right)$
6	3	$1228.78$	$\left(1742.74\right)$	80	$0.53$	$\left(2.55\right)$
7	0	-		86	$0.74$	$\left(4.92\right)$
8	0	-		94	$0.86$	$\left(4.39\right)$
9	0	-		96	$0.90$	$\left(3.39\right)$
10	0	-		94	$1.47$	$\left(8.69\right)$
11	0	-		97	$2.45$	$\left(13.97\right)$
12	0	-		99	$4.86$	$\left(22.97\right)$
13	0	-		98	$9.64$	$\left(56.71\right)$
14	0	-		100	$24.31$	$\left(84.90\right)$
15	0	-		100	$60.04$	$\left(236.27\right)$
16	0	-		100	$245.58$	$\left(1291.61\right)$

Equations

$$\text{Find }x\in \bigcap _{i=1}^N{C}_i.$$

$$P_C(x):=\left\{z\in C:\Vert x-z\Vert =\underset{c\in C}{\inf }\Vert x-c\Vert \right\},$$

$$p=P_C(x)⟺p\in C\text{ and }⟨c-p,x-p⟩\le 0,\forall c\in C.$$

$$x_{n+1}=T_{A,B}(x_n)\text{for }n=0,1,2,\dots .$$

$$\overline{x}=2P_B(2P_A(\overline{x})-\overline{x})-(2P_A(\overline{x})+\overline{x})⟺P_A(\overline{x})=P_B(2P_A(\overline{x})-\overline{x}),$$

$$\Vert P_A(x)-P_A(y){\Vert }^2+\Vert (\mathrm{Id}-P_A)(x)-(\mathrm{Id}-P_A)(y){\Vert }^2\le \Vert x-y{\Vert }^2,\forall x,y\in \mathcal{H},$$

$$\Vert R_A(x)-R_A(y)\Vert \le \Vert x-y\Vert ,\forall x,y\in \mathcal{H},$$

$${\mathcal{H}}^N:=\mathcal{H}\times \stackrel{(N)}{\cdots }\times \mathcal{H}=\{\mathbf{x}=(x_1,\dots ,x_N):x_i\in \mathcal{H},i=1,\dots ,N\}.$$

$$C:=C_1\times C_2,\times \cdots \times C_N\text{and}D:=\left\{(x,x,\dots ,x)\in {\mathcal{H}}^N:x\in \mathcal{H}\right\}.$$

$$x\in \bigcap _{i=1}^N{C}_i⟺(x,x,\dots ,x)\in C\cap D.$$

$$P_D(\mathbf{x})=\left(\frac{1}{N}\sum _{i=1}^N{x}_i,\dots ,\frac{1}{N}\sum _{i=1}^N{x}_i\right),$$

$$P_C(\mathbf{x})=P_{C_1}(x_1)\times P_{C_2}(x_2)\times \cdots \times P_{C_N}(x_N).$$

$$0=⟨\mathbf{d},\mathbf{x}-\mathbf{p}⟩=\sum _{i=1}^N⟨d,x_i-p⟩=⟨d,\sum _{i=1}^N{x}_i-Np⟩.$$

$$\Vert \mathbf{x}-\mathbf{c}{\Vert }^2=\sum _{i=1}^N\Vert x_i-c_i{\Vert }^2\ge \sum _{i=1}^N\Vert x_i-p_i{\Vert }^2=\Vert \mathbf{x}-\mathbf{p}{\Vert }^2.$$

$$\sum _{i=1}^N\Vert x_i-p_i{\Vert }^2=\Vert \mathbf{x}-\mathbf{p}{\Vert }^2\le \Vert \mathbf{x}-\mathbf{q}{\Vert }^2=\sum _{i=1}^N\Vert x_i-q_i{\Vert }^2=\sum _{\genfrac{}{}{0ex}{}{i=1}{i\ne j}}^N\Vert x_i-p_i{\Vert }^2+\Vert x_j-c_j{\Vert }^2,$$

$$x_{n+1}=T_{C_N,C_1}{T}_{C_{N-1},C_N}\cdots T_{C_2,C_3}{T}_{C_1,C_2}(x_n);$$

$$C_1$$

$$C_2$$

$$C_3$$

$$C_4$$

$$P_{\mathcal{C}}(x)=[\mathcal{C}{]}_x,$$

$$B_{i,j,k}:=\{\begin{array}{ll}1& \text{if }{A}_{i,j}=k,\\ 0& \text{otherwise.}\end{array}$$

$$C_1$$

$$C_2$$

$$C_3$$

$$C_4$$

$$C_5=\left\{B\in {\mathbb{R}}^{9\times 9\times 9}:B_{i,j,:}\in \mathcal{B}\text{ for each }i,j\in I\right\}.$$

$$C_1$$

$$C_2$$

$$C_3$$

$$C_4$$

$$C_5$$

$$P_S(x)=x+\frac{1}{m}\left(c-\sum _{i=1}^m{x}_i\right)e,\text{where }e=[1,1,\dots ,1{]}^T.$$

$$P_{C_1}(A)$$

$$P_{C_2}(A)$$

$$P_{C_3}(A)$$

Full text

Keywords: Douglas–Rachford algorithm, magic square, feasibility problem, nonconvex constraints

AMS Subject classifications: 47J25, 47N10, 05B15, 97A20

1 Introduction

A magic square of order $n$ is an $n\times n$ matrix whose elements are distinct positive integers $1,2,\ldots,n^{2}$ arranged in such a way that the sum of the $n$ numbers in any horizontal, vertical, or main diagonal line is always the same number, known as the magic constant. Clearly, the magic constant is equal to $n(n^{2}+1)/2$, as it is the result of dividing by $n$ the sum of $1,2,\ldots,n^{2}$. Hence, for magic squares of order $3$, $4$ and $5$, their respective magic constants are $15$, $34$ and $65$. In Fig. 1 we show the unique magic square of order $3$, except for rotations and reflections.

Magic squares have a long history and were frequently associated with mystical and supernatural properties. It is believed that the magic square of order $3$ shown in Fig. 1, commonly referred by its Chinese name luoshu, was known to Chinese mathematicians long before 500 BCE. According to [26], the very first textual reference to the luoshu appears to be in the writings of Zhuang Zi (369-286 B.C.E.), one of the founders of Daoism. The name luoshu means “Luo River writing”, and refers to the following legend:

There was a huge flood in the ancient China. While Sage King Yu stood on the banks of the Luo River trying to channel the water out to the sea, a turtle with a curious pattern of dots arranged on its shell emerged from the river. This pattern was the magic square of order 3. Thereafter people were able to use this pattern in a certain way to control the river and protect themselves from floods.

It is possible to construct magic squares of any order, except for $n=2$: if the following grid was a magic square,

{TAB}

(e,0.6cm,0.6cm)|c|c||c|c| a & b

c d

then one should have $a+b=a+c$, which implies $b=c$.

Probably, the most well-known magic square of order $4$ is the one immortalized by the German artist Albrecht Dürer in his engraving Melencolia I, see the detail in Fig. 2. Observe how the date of the engraving, 1514, is shown in the two middle cells of the bottom row. This magic square has the additional property that the sums in any of the four quadrants, as well as the sum of the middle four numbers, are all $34$ (which is called a gnomon magic square).

Dürer became very interested in mathematics and devoted his life to the study of linear perspective and proportion. He claimed:

The new art must be based upon science – in particular, upon mathematics, as the most exact, logical, and graphically constructive of the sciences.

Dürer’s magic square previously appeared in the manuscript De viribus quantitatis by the Franciscan friar Luca Pacioli, which includes a short chapter on magic squares dating from about 1501–1503. In that chapter, Pacioli includes one example each of magic square of orders $3$ to $9$, and the one of order $4$ is identical to the one in Melancolia I. Note that the chances for coincidence are very low, as there are $880$ magic squares of order $4$ (multiplied by their rotations and reflections). It is very plausible that Dürer learnt this magic square from Pacioli, during his second visit to Italy between 1505–1507. For more details, see [26, 27].

There are many ways of constructing magic squares. In this work, based on [24], we present a novel approach by formulating their construction as a nonconvex feasibility problem which can be tackled by the Douglas–Rachford algorithm. Of course, this only an expository problem, as there are much more efficient strategies for constructing magic squares. We review the main characteristics of this algorithm in Section 2. Two different formulations are presented in Section 3. The numerical experiments reported in Section 4 demonstrate that the Douglas–Rachford algorithm can be used as a heuristic for finding them. We finish with some conclusions in Section 5.

2 Projection algorithms

Given a finite family $C_{1},C_{2},\ldots,C_{N}$ of subsets of a Hilbert space $\mathcal{H}$, a feasibility problem consists in finding a point in their intersection; that is,

[TABLE]

Thanks to the generality of (2.1), plenty of problems can be modeled as feasibility problems. In many cases, finding a point in the intersection of the sets is a challenging task, but computing the projection of points into each of the sets is easy. In such scenarios, projection methods become very useful. These are algorithms that use the projection of points onto the sets to define the iterates. Recall that the projection of a point $x$ onto a closed set $C\subset\mathcal{H}$ is defined by

[TABLE]

where $\|\cdot\|$ denotes the norm induced by the inner product $\langle\cdot,\cdot\rangle$ of $\mathcal{H}$. Observe that, in general, the projection operator $P_{C}:\mathcal{H}\rightrightarrows C$ is a set-valued mapping, see Fig. 3. Nonetheless, if the set $C$ is nonempty, closed and convex, the projection is uniquely determined as follows (see, e.g., [8, Theorem 3.16]):

[TABLE]

The most intuitive and well-known projection method is the von Neumann’s alternating projection algorithm, which iterates by cyclically projecting onto each of the sets, see Fig. 4.

In this work, we will focus on another popular projection method: the Douglas–Rachford algorithm. Although the scheme has its origins in the work of J. Douglas and H.H. Rachford [14], where it was proposed for solving a system of linear equations arising in heat conduction problems, it was P.L. Lions and B. Mercier [19] who really deserve the credit for the algorithm. In their work, not only they showed how to successfully generalize the algorithm for solving convex feasibility problems, but they actually provided a splitting algorithm for finding a zero in the sum of two maximally monotone operators.

The Douglas–Rachford algorithm employs the reflection operators when is applied for solving feasibility problems. Recall that the reflection of a point $x\in\mathcal{H}$ onto a closed set $C$ is defined as $R_{C}(x):=2P_{C}(x)-x$. The iterates of the algorithm for solving feasibility problems involving two sets are obtained by computing an average between the current point and two consecutive reflections onto the sets, see Fig. 5. For this reason, the scheme is also referred as the averaged alternating reflections method.

The Douglas–Rachford algorithm has recently gained great attention, in part thanks to its good performance as a heuristic in nonconvex settings (where the convergence of the algorithm is not guaranteed). In the remainder of this section we will review its main properties in the convex setting, we will show how the algorithm can be applied to any number of sets thanks to a product space reformulation of the feasibility problem, and we will review an application in the nonconvex setting for solving Sudoku puzzles.

2.1 The Douglas–Rachford algorithm for convex sets

When the Douglas–Rachford algorithm is applied to two closed and convex sets $A,B\subset\mathcal{H}$, it can be viewed as a fixed point iterative method defined by the Douglas–Rachford operator $T_{A,B}:=\frac{1}{2}(\operatorname{Id}+R_{B}R_{A})$, where $\operatorname{Id}$ denotes the identity mapping:

[TABLE]

Note that the iterates are uniquely determined because the reflection operators onto convex sets are single-valued, as the projection operators are. Further, observe that $\overline{x}\in\mathcal{H}$ is a fixed point of the operator $T_{A,B}$ if and only if

[TABLE]

which implies that $P_{A}(\overline{x})\in A\cap B$. Therefore, the projection onto the set $A$ of any fixed point of the Douglas–Rachford operator solves the feasibility problem.

Using a well-known theorem of Opial [20, Theorem 1], the sequence generated by the algorithm can be proved to be weakly convergent to a fixed point of the Douglas–Rachford operator, whenever $A\cap B\neq\emptyset$. This is a consequence of the firm nonexpansivity of the projection operator,

[TABLE]

which implies that the reflection operator is nonexpansive (or Lipschitz continuous with constant $1$), i.e.,

[TABLE]

and therefore, $T_{A,B}$ is firmly nonexpansive [17, Theorem 12.1].

The next result summarizes the main properties of the Douglas–Rachford algorithm in the convex setting.

Theorem 2.1**.**

Let $A,B\in\mathcal{H}$ be closed and convex sets. Given any $x_{0}\in\mathcal{H}$, define $x_{n+1}=T_{A,B}(x_{n})$, for every $n\geq 0$. Then, the following holds:

• (i)

If $A\cap B\neq\emptyset$, then $\{x_{n}\}$ is weakly convergent to a point $x^{\star}$ and $\{P_{A}(x_{n})\}$ is weakly convergent to $P_{A}(x^{\star})\in A\cap B$,

• (ii)

If $A\cap B=\emptyset$, then $\|x_{n}\|\rightarrow\infty$.

Proof. (i) For the first part, see [9, Theorem 3.13 and Corollary 3.9]. Note that the weak convergence of the shadow sequence $\{P_{A}(x_{n})\}$ cannot be derived from the weak convergence of $\{x_{n}\}$, since $P_{A}$ may not be weakly continuous (see [8, Example 4.20]). The weak convergence of the shadow sequence was proved by Svaiter in [25, Theorem 1]. (ii) See [7, Corollary 2.2].

In contrast with the alternating projection algorithm, whose scheme can be easily generalized to more than two sets by cyclically projecting onto the sets, the same cannot be done with the Douglas–Rachford algorithm (see [1, Example 2.1] for a simple example involving three lines in $\mathbb{R}^{2}$). Fortunately, the product space formulation, due to Pierra [22], permits to reformulate any feasibility problem as another one involving only two sets. This is the subject of the next subsection.

2.2 Product space formulation

Consider the feasibility problem (2.1) given by $N$ sets. Let

[TABLE]

It can be easily to checked that the space $\mathcal{H}^{N}$ endowed with the inner product $\langle\mathbf{x},\mathbf{y}\rangle:=\sum_{i=1}^{N}\langle x_{i},y_{i}\rangle$ is a Hilbert space. Let us consider two sets $C,D\subset\mathcal{H}^{N}$ defined by

[TABLE]

The set $D$, which is sometimes called the diagonal, is always a closed subspace. Observe that the feasibility problem (2.1) is equivalent to the one involving the sets $C$ and $D$, since

[TABLE]

Not only (2.1) can be recast as a feasibility problem involving two sets, but also the projection onto each of these sets has a closed form in terms of $P_{C_{1}},P_{C_{2}},\ldots,P_{C_{N}}$, as shown next. For completeness, we include a proof based on [2, Proposición 3.1]. This result extends Pierra’s Lemma 1.1 in [22] to possibly nonconvex sets.

Proposition 2.1**.**

For any $\mathbf{x}=(x_{1},x_{2},\ldots,x_{N})\in\mathcal{H}^{N}$, one has

[TABLE]

and if $P_{C_{1}}(x_{1}),P_{C_{2}}(x_{2})\ldots,P_{C_{N}}(x_{N})$ are nonempty,

[TABLE]

Proof. Let $\mathbf{p}=(p,\ldots,p)\in D$ be the projection of $\mathbf{x}$ onto the set $D$. For any $d\in\mathcal{H}$, one has $\mathbf{d}:=(d,\ldots,d)\in D$. Then, using the characterization (2.2) applied to $\mathbf{p}+\mathbf{d}$ and $\mathbf{p}-\mathbf{d}$, we deduce

[TABLE]

As the latter holds for any $d\in\mathcal{H}$, it must be $p=\frac{1}{N}\sum_{i=1}^{N}x_{i}$, as claimed.

To prove the formula for the projection onto the set $C$, pick first any point $\mathbf{p}=(p_{1},\ldots,p_{N})\in P_{C_{1}}(x_{1})\times\cdots\times P_{C_{N}}(x_{N})\subseteq C$. Then, for any $\mathbf{c}=(c_{1},\ldots,c_{N})\in C$, one has

[TABLE]

Therefore, $P_{C}(\mathbf{x})\supseteq\prod_{i=1}^{N}P_{C_{i}}(x_{i})$.

Conversely, let $\mathbf{p}=(p_{1},\ldots,p_{N})\in P_{C}(\mathbf{x})$. Pick any $j\in\{1,\ldots,N\}$ and choose any $c_{j}\in C_{j}$. Let $q_{i}:=p_{i}$ if $i\neq j$ and let $q_{j}:=c_{j}$, and define $\mathbf{q}:=(q_{1},\ldots,q_{N})\in C$. Then,

[TABLE]

whence, $\|x_{j}-p_{j}\|\leq\|x_{j}-c_{j}\|$. Since $c_{j}$ was arbitrarily chosen in $C_{j}$, this implies that $p_{j}\in P_{C_{j}}(x_{j})$, and concludes the proof.

Pierra’s product space formulation works well in practice whenever the number of sets $N$ is small. When $N$ is large, the computational cost of calculating the iterates in the large space $\mathcal{H}^{N}$ instead of $\mathcal{H}$ becomes prohibitive. In these cases, the cyclic Douglas–Rachford algorithm with $r$-sets-Douglas–Rachford operators, recently introduced in [6], would be a better choice. The latter is a generalization of the cyclic Douglas–Rachford method of Borwein–Tam [13], whose iterates are defined by

[TABLE]

that is, a cyclic composition of Douglas–Rachford operators applied to the sets.

2.3 The Douglas–Rachford algorithm for nonconvex sets

Theorem 2.1 only guarantees the global convergence of the Douglas–Rachford algorithm for convex sets. In spite of this, the algorithm has been successfully applied as a heuristic for solving many nonconvex problems, particularly those of combinatorial type (see, e.g., [1, 2, 4, 10, 15]). In most applications in the nonconvex setting, the constraint sets satisfy some type of regularity property and local convergence can be proved [11, 18, 21]. Nonetheless, the results on global behavior are limited to very special sets [3, 12] and cannot explain the good performance of the algorithm in the nonconvex setting.

It is clear that the performance of an algorithm heavily depends on how the problem is formulated. When the Douglas–Rachford algorithm is applied to a nonconvex problem, finding an “adequate” formulation can be crucial for its success as a heuristic. Indeed, the formulation of the problem often determines whether or not the algorithm can successfully solve the problem at hand always, frequently or never.

Let us briefly present an interesting example exhibiting this behavior: Sudoku puzzles. Their solution with the Douglas–Rachford algorithm was first proposed in [15], and subsequently analyzed in [23] and [2]. A Sudoku is a $9\times 9$ grid, divided in nine $3\times 3$ subgrids, with some entries already prefilled. The objective is to fill the remaining entries in such a way that each row, column and subgrid contains the digits from $1$ to $9$ exactly once, see Fig. 6.

Then, Sudokus are nothing else than matrix completion problems, which can be modeled as feasibility problems. There are different possibilities for choosing the constraint sets $C_{1},\ldots,C_{N}$ in such a way that $\bigcap_{i=1}^{N}C_{i}$ coincides with the (unique) solution to the Sudoku. Nonetheless, in order to apply the Douglas–Rachford algorithm, the sets must be chosen in such a way that the projections can be efficiently computed, ideally having a closed form. In [2, Sección 6], two ways of modeling Sudokus are analyzed: as integer and binary feasibility problems. We briefly describe them in what follows.

Sudoku modeled as an integer program

It is obvious how to model a Sudoku as an integer feasibility problem. Let $S$ be the partially filled $9\times 9$ matrix representing the unsolved Sudoku, let $J\subset\{1,2,\ldots,9\}^{2}$ be the set of indices for which $S$ is filled, and let $A_{i,j}$ denote the $(i,j)$-th entry of the matrix $A$. If we denote by $\mathcal{C}$ the set of vectors which are permutations of $1,2,\ldots,9$, then a matrix $A\in\mathbb{R}^{9\times 9}$ is a solution to the Sudoku if and only if $A\in C_{1}\cap C_{2}\cap C_{3}\cap C_{4}$, where

[TABLE]

The projection onto $C_{4}$ can be easily computed componentwise, while the projections onto $C_{1}$, $C_{2}$ and $C_{3}$ are determined by the next result, see [2] for more details.

Proposition 2.2**.**

Denote by $\mathcal{C}\in\mathbb{R}^{m}$ the set of vectors whose entries are all permutations of $c_{1},c_{2},\ldots,c_{m}\in\mathbb{R}$. Then, for any $x\in\mathbb{R}^{m}$, one has

[TABLE]

where $[\mathcal{C}]_{x}$ denotes the set of vectors in $\mathcal{C}$ such that $y\in[\mathcal{C}]_{x}$ if the $i$th largest entry of $y$ has the same index in $y$ as the $i$th largest entry of $x$.

Proof. See [2, Proposition 5.1].

Sudoku modeled as a zero-one program

To model Sudoku puzzles as binary programs, we reformulate a matrix $A\in\mathbb{R}^{9\times 9}$ as $B\in\mathbb{R}^{9\times 9\times 9}$, where

[TABLE]

This reformulation transforms the entries into a 3-dimensional zero-one array, which can be thought as a cube, see Fig. 7.

Let us denote by $S^{\prime}$ the $9\times 9\times 9$ zero-one array corresponding to the unsolved Sudoku $S$, let $I:=\{1,2,\ldots,9\}$, let $J^{\prime}\subseteq I^{3}$ be the set of indices for which $S^{\prime}$ is filled, and let $\mathcal{B}$ denote the $9$-dimensional standard basis. The four constraints of the integer program above become now

[TABLE]

where $\operatorname{vec}A$ represents the vectorization of a matrix $A$ by columns.

Further, since each cell in the Sudoku can only have one number assigned, we must add the additional constraint

[TABLE]

These constraints are represented in Fig. 8.

Then, $B$ completes $S^{\prime}$ (and thus solves the Sudoku) if and only if $B\in C_{1}\cap C_{2}\cap C_{3}\cap C_{4}\cap C_{5}$. Again, the projection onto $C_{4}$ can be easily computed componentwise, while the projections onto $C_{1}$, $C_{2}$, $C_{3}$ and $C_{5}$ are determined by Proposition 2.2.

Performance of the Douglas–Rachford algorithm on Sudokus

As observed in [2], the Douglas–Rachford algorithm is totally ineffective for solving the integer formulation. On the other hand, the algorithm is very successful when it is applied to the binary formulation, being able to solve nearly all the instances (without a restart) in all the Sudoku libraries tested in [2]. Even in the one where it was most unsuccessful, the top95111top95: http://magictour.free.fr/top95 library, it solved 87% of the instances.

Another possibility for modeling Sudoku as feasibility problems was tested in [4, 5], based on the fact that an unsolved Sudoku can be viewed as a graph precoloring problem. Surprisingly, the rank formulation proposed in [5] had a success rate of 100% in the top95 library, and no Sudoku has been found so far for which Douglas–Rachford fails to find its solution for any starting point.

Therefore, it is clear that the choice of the formulation has a big impact in the success rate of the Douglas–Rachford algorithm when it is applied in the nonconvex setting.

3 Finding magic squares as a feasibility problem

In this next section, we propose two different formulations for finding magic squares as a feasibility problem which are inspired by the binary and integer formulations of Sudokus. To formulate the search, the first thing we need to do is to choose some sets whose intersection includes all the properties that characterize a magic square of order $n$: the sum of each row, column and main diagonals must be equal to the magic constant $c:=\frac{n(n^{2}+1)}{2}$, and it must contain all the numbers between $1$ and $n^{2}$. Note that if our purpose was to complete a prefilled magic square, we would need to add an additional constraint that fixes the prefilled numbers.

The most obvious way to look for a magic square is to use an integer formulation. This is explained in the next subsection, where we also deduce the expression for the projection onto each of the constraint sets. Afterwards, we show how to formulate the search for magic squares as a zero-one program, in line with what we have seen for Sudokus in Section 2.3, although in a slightly different way.

3.1 Magic squares modeled as an integer program

If we denote by $\mathcal{P}$ the set of permutations of $1,2\ldots,n^{2}$ and $I:=\{1,2,\ldots,n\}$, then $A\in\mathbb{R}^{n\times n}$ is a magic square if and only if $A\in\bigcap_{i=1}^{5}C_{i}$, where

[TABLE]

The first four sets are clearly convex, and their projection operators have a closed form determined by the next result.

Proposition 3.1**.**

Consider $S=\left\{x\in\mathbb{R}^{m}:\sum_{i=1}^{m}x_{i}=c\right\}$. For any $x\in\mathbb{R}^{m}$,

[TABLE]

Proof. This follows from the standard formula for the orthogonal projection onto a hyperplane, since $S=\left\{x\in\mathbb{R}^{m}:\langle x,e\rangle=c\right\}$; see, e.g., [8, Example 3.23].

Thereby, the projections onto each of these sets are given by

[TABLE]

where $I_{n}$ denotes the identity matrix in $\mathbb{R}^{n\times n}$.

Finally, $C_{5}$ is a nonconvex set containing all those matrices whose entries are permutations of $1,2,\ldots,n^{2}$, so its projection operator is determined by Proposition 2.2. Therefore, given an $n\times n$ matrix, to compute its projection onto $C_{5}$, we just need to sort the elements of the matrix in ascending order and place the number $1$ in the cell that contains the smallest number, the number $2$ in the cell that contains the next one, and so on. If there are two equal elements, the projection of the matrix is not unique (in our experiments in Section 4, we chose the element that appeared first when the matrix was read by rows).

Completing a partially filled magic square

If $M$ is a partially complete matrix representing an incomplete magic square, denote by $J\subseteq I^{2}$ the set of indices for which $M$ is filled. To find a completion of the magic square, we need to add the constraint

[TABLE]

to the set of constraints that we already have. Then, $A$ completes $M$ if, and only if,

[TABLE]

The projection onto $C_{6}$ is given componentwise by

[TABLE]

for each $(i,j)\in I^{2}$.

3.2 Magic squares modeled as a binary program

To model the search for a magic square as a binary feasibility problem, we reformulate a matrix $A\in\mathbb{R}^{n\times n}$ as $B\in\mathbb{R}^{n\times n\times n^{2}}$, where

[TABLE]

In this way, we transform the entries of the matrix into a $3$-dimensional zero-one array, and each number in the magic square can be thought as a pillar made of small cubes, see Fig. 9.

The constraints from the previous subsection become (see Fig. 10):

[TABLE]

Constraint $C_{5}$ of the integer formulation is now written for the binary formulation as the intersection of two sets:

[TABLE]

where $\mathcal{B}:=\left\{[1,1,\ldots,1,1],[1,1,\ldots,1,0],\ldots,[1,0,\ldots,0,0]\right\}$ is a base of vectors of $\mathbb{R}^{n^{2}}$. On the one hand, constraint $C_{5}$ guarantees that the first floor of $B$ is all filled with ones, the second must have $n^{2}-1$ ones, and so on until we get to the top floor, which only contains a one. On the other hand, constraint $C_{6}$ guarantees that the matrix is formed by pillars of ones standing on the floor, so if there is a one in an entry, all the elements below must be ones as well. Then, $B\in\mathbb{R}^{n\times n\times n^{2}}$ is a magic square if and only if $B\in\bigcap_{i=1}^{6}C_{i}.$

The projections onto the first four sets are given by Proposition 3.1:

[TABLE]

where $e=(1,\ldots,1)\in\mathbb{R}^{n^{2}}$ and

[TABLE]

The projection onto $C_{5}$ assigns a one to all the elements of the first floor, a zero to the smallest element and a one to the rest in the second floor, a zero to the two smallest elements and a one to the rest in the third floor, and so on until the top floor, which only contains a one. The projection onto $C_{6}$ can be computed componentwise as follows: For each $i,j\in I$, we have

[TABLE]

Remark 3.1**.**

The formulation of the problem of completing a partially filled magic square is analogous to the integer formulation, so we omit the details.

4 Numerical experiments

In this section we run a numerical experiment to test the performance of the Douglas–Rachford algorithm for finding magic squares of order $n=3,4,\ldots,12$, on which we compare the integer and binary formulations discussed in Section 3. All our codes were written in Python 2.7 and the tests were run on an Intel Core i7-4770 CPU 3.40GHz with 32GB RAM, under Windows 10 (64-bit).

For each order and each formulation, the algorithm was run on the corresponding product space (see Section 2.2) for $100$ random starting points. For the integer formulation, we took a starting point of the form $x_{0}:=(y,y,y,y,y)\in D$ for some $y\in\mathbb{R}^{n\times n}$ randomly chosen with entries in $]0,1[$, while for the binary formulation, we choose $x_{0}:=(y,y,y,y,y,y)\in D$ for some random $y\in\mathbb{R}^{n\times n\times n^{2}}$ with entries in $]0,1[$. The iterates were then defined by

[TABLE]

The Douglas–Rachford algorithm was stopped either after a maximum time of $1800$ seconds ($30$ minutes) or when

[TABLE]

where $\operatorname{round}(\cdot)$ gives the nearest integer componentwise.

Table 1 summarizes the success rate and running time for each formulation. We observe how the binary formulation was not able to solve any instance for orders larger than $6$. On the other hand, for orders $4$ and $5$, it was more successful than the integer formulation, which surprisingly becomes more successful for larger orders than for smaller ones. The algorithm was much faster with the integer formulation.

To better understand the behavior of the Douglas–Rachford algorithm for both formulations, we show in Fig. 11 the cumulative frequency over time. Note that we have used a logarithmic scale in the horizontal axis, to show how the running time of the algorithm exponentially grows with the order $n$. Apparently, the reason why the algorithm failed in the binary formulation for $n$ larger than $6$ is that it did not have enough time to converge. We show in Fig. 12 a magic square of orders between $3$ and $12$ found with the integer formulation.

5 Conclusions

In this article, we have shown that the Douglas–Rachford algorithm can be used as a successful heuristic for constructing magic squares, despite the nonconvexity of the problem. To this aim, we have proposed two different formulations as a feasibility problem, one integer and another one binary. Our numerical tests demonstrate that the integer formulation is more effective for finding magic squares because it requires considerably much less time on average for finding a solution. Surprisingly, this behavior was the opposite with the similar binary and integer formulations for solving Sudoku puzzles, where the integer formulation is completely unsuccessful. Interestingly, the binary was more successful than the integer formulation for orders $4$ and $5$. For higher orders, the binary formulation was not able to solve any instance, possibly due to lack of time.

Acknowledgements

The first author was supported by MINECO of Spain and ERDF of EU, as part of the Ramón y Cajal program (RYC-2013-13327) and the Grant MTM2014-59179-C2-1-P.

About the authors

Francisco J. Aragón Artacho earned his Ph.D. in 2007 from the University of Murcia (Spain). He subsequently held a postdoctoral position “Juan de la Cierva” at the University of Alicante. In 2011 he moved to the University of Newcastle (Australia), where he became a Research Associate at the Research Centre for Computer-Assisted Research Mathematics and its Applications. Two years later, he moved to Luxembourg where he became a Marie Curie Fellow (AFR Postdoc) at the Luxembourg Centre for Systems Biomedicine. In November 2014 he started his current position at the University of Alicante as a Ramón y Cajal fellow. His main research interests include set-valued and variational analysis, optimization, convex analysis and projection algorithms.

Paula Segura Martínez graduated from the University of Alicante with a B.S. in Mathematics in 2017. She received her M.S. in Mathematics from the University of Murcia in 2018 and now she is a Ph.D. student at the University of Valencia. Her research interests are in the fields of combinatorial optimization, location and arc routing problems.