OPTIMIZATION-BASED APPROACH TO TILING OF FINITE AREAS WITH ARBITRARY SETS OF WANG TILES

Wang tiles proved to be a convenient tool for the design of aperiodic tilings in computer graphics and in materials engineering. While there are several algorithms for generation of finite-sized tilings, they exploit the specific structure of individual tile sets, which prevents their general usage. In this contribution, we reformulate the NP-complete tiling generation problem as a binary linear program, together with its linear and semidefinite relaxations suitable for the branch and bound method. Finally, we assess the performance of the established formulations on generations of several aperiodic tilings reported in the literature, and conclude that the linear relaxation is better suited for the problem.


Introduction
Wang tiles, squares with colored edges, were invented by and named in honor of Hao Wang, originally serving as a tool for studying the ∀∃∀ decidability problem of the predicate calculus [1].Wang showed that the decidability problem is equivalent to the domino problem: assume a set of non-rotatable unit-sized (Wang) tiles with edges colored according to the ∀∃∀ problem.If it is feasible to cover the infinite plane by translating copies of the tiles, such that all their vertices lie at the integer lattice points of the plane and the adjoining edges share the same color, the problem is said to be solvable; unsolvable otherwise.For an illustration of a sample 2 × 3 Wang tiling, see Fig. 1.
In [2], Wang made a fundamental conjecture stating that the domino problem is solvable if and only if the tiling was periodic, i.e., if there existed a rectangular region of the tiling with identically colored horizontal, and vertical edges, respectively.
A year later, Kahr reduced the Turing machine halting problem [3,4] into the origin-constrained domino problem [5], which implies the domino problem is also undecidable.This can be illustrated by introducing a Turing machine for each tile set, halting only if the domino problem is unsolvable.However, there is an infinite number of such tile sets, making the domino problem undecidable as well, and forbidding existence of a general finite algorithm for the generation of infinite tilings.

Aperiodic Tile Sets
Undecidability of the domino problem was also proved by Wang's student Berger, who used the principle of expanding squares for developing the sets of 20, 426 [6] and 104 tiles [7] that admit only aperiodic tilings of the infinite plane, contrary to the Wang fundamental conjecture.Aiming to simplify the Berger's proof, smaller aperiodic sets of Wang tiles have been constructed.In 1966, Läuchli sent to professor Wang an aperiodic set of 40 tiles over 16 colors, but it remained unpublished until 1975 [8].Meanwhile, unaware of Läuchli's result, Knuth simplified Berger's set to 92 tiles over 26 colors [9]; and Robinson developed sets of 104 and 52 tiles in 1967 [10], of 56 tiles over 12 colors in 1971 [11], and marked existence of a set of 35 tiles [11].
In 1973, Penrose developed a new approach based on kites and darts tiling, leading to a set of 34 tiles.Robinson, being in contact with Penrose, modified the Penrose's approach to reach a reduced set of 32 tiles over 16 colors [12].Using the same technique together with Penrose rhombs tiling, Grünbaum obtained a set of 24 tiles over 9 colors in 1987 [12].
Another two tile sets were discovered due to Ammann.In 1978, he used Ammann bars to reach 16 tiles over 6 colors [13].Building on the Ammann's A2 tiling, see, e.g., [12], Robinson obtained a set of 24 tiles over 24 colors in 1977.
Subsequent size reduction of the smallest aperiodic set occurred in 1996, as Kari developed a new method based on Mealy machines multiplying Beatty sequences, and presented a set of 14 tiles over 6 colors [14].Čulík, using the same approach, reduced the set even further to 13 tiles over 5 colors [15].
The search for the minimal aperiodic set is concluded by Jeandel, who used a brute-force enumeration approach to find aperiodic sets of 11 tiles over 4 and 5 colors; and proved both that there does not exist an aperiodic set with 10 or fewer tiles and with less than 4 colors [16].
In addition to the classical Wang tiles, in 2006, Lagae introduced a subset of the Wang tiles, the corner tiles, with the connectivity information stored in the colored corners instead of the edges [17].In the same year, he constructed an aperiodic set of 44 corner tiles over 6 colors [18].The set was further simplified by Nurmi to 30 corner tiles over 6 colors [19].Note that because the corner tiles can be straightforwardly converted to the edge-based Wang tiles [17], our approach naturally extends to the corner tiles, see Section 4 for a specific example.

Selected Applications
Due to the ability of some tile sets to generate aperiodic patterns, Wang tiles received a broad interest among disciplines.Also, their original significance for automated theorem proving [2] was supplemented by proofs in cellular automata theory [20].
In 2003, Cohen et al. introduced Wang tiles to the computer graphics community [21].Since then, Wang tiles have been successfully used, e.g., for efficient synthesis of stochastic textures or for generation of Poisson disk distributions.
Molecular realization of Wang tiles is due to Winfree, who introduced a self-assembly of biological nanostructures into aperiodic crystals [22].
Novák was probably the first one who used Wang tiles in materials engineering for efficient compression of microstructures [23]; Doškář for their reconstruction [24].Doškář further generated large stochastic samples of the Alporas® foam and its finiteelement representations in [25].In [26], Tyburec used Wang tiles to describe modular assembly of structures, whereas both the topology and arrangement of modules were subjects of optimization.

Tiling Generation Algorithms
Existing tiling generation algorithms are designed specifically to the tile sets they handle.In computer graphics, for example, it is essential to generate visually appealing patterns quickly, which is best achieved with stochastic tile sets.Tiling a finite area has then a O(n) complexity, adopting, e.g., the stochastic tiling [21], or the hash-based direct stochastic tiling [27] algorithms.The latter algorithm allows straightforward definition of boundary conditions.
In the case of mathematical logic, it is crucial to prove that an investigated tile set can tile the whole infinite plane, and does so only aperiodically.As the problem is undecidable, no general algorithm exists, and thus all the algorithms need to exploit the specific structure of each individual (family of) tile sets.Although the algorithms are fast, they remain tile set dependent.All tilings of finite-sized areas are restricted to be subsets of the infinite plane, which is not required in the finite domain, and it is almost impossible to introduce boundary conditions.On top of that, there provably exist aperiodic tile sets that admit only non-recursive tiling [28,29], forbidding design of any problem-specific algorithm to tile the infinite plane.
To the authors' knowledge, the only method that has been used for tiling of finite-sized areas by arbitrary tile sets is the backtracking algorithm, see [30].Although the algorithm is general, it is generally inefficient, as it commonly creates impossible assemblies too early.Moreover, distant boundary conditions make the problem more difficult to solve.

Contributions
This paper aims to overcome the shortcomings of the methods outlined in the previous section, and to develop an approach that handles tiling of finitesized areas using arbitrary tile sets, together with a straightforward approach to define edge-or tilebased boundary conditions.
Our exposition is structured as follows.In Section 2, we develop a binary linear programming representation of the tiling generation problem.The formulation is relaxed, and its linear and semidefinite approximations are presented.In order to show generality of the formulations, we introduce their extensions to solve the tile-packing problem, to prove that a tile set can not tile the plane, and to prove that a tile set admits periodic tiling.The relaxations are finally employed in a branch and bound method in Section 3 and their performance is assessed in Section 4.

Valid Tiling
Consider a finite rectangular area A of the size n t,h × n t,w , with n t,h and n t,w denoting its height and width, respectively.Placing Wang tiles (of a unit size) from the tile set T , such that their vertices lie at the integer lattice points of A, we obtain a tiling with n t,h rows and n t,w columns of tiles.Each tile k ∈ {1, .. , n t } from the tile set T , with n t denoting the number of tiles, is described by the 4-tuple (n k , w k , s k , e k ), assigning color codes C to the north, west, south, and east edge of the tile, respectively.
In any valid tiling, for any two adjacent tiles k and , with the tile being placed in the east of k, it holds that e k = w . ( Similarly, for any two vertically adjacent tiles k and m, with the tile m being placed in the south of k, it stands that Both the cases are illustrated in Fig. 2.

Binary Linear Programming Formulation
Let x ∈ {0, 1} n t,h •nt,w•nt denote a binary vector describing the placement of individual tiles within A.
The vector consists of all the x i,j,k variables, where i ∈ {1, .. , n t,h } and j ∈ {1, .. , n t,w } are the row and column iterators, respectively.Let us also define Because each position is occupied by a single tile, we have The compatibility constraint (1) can be written in terms of x after realizing that the number of tiles at (i, j) that have east edge colored by c ∈ C has to be equal to the number of tiles at (i, j + 1) with west edge colored also by c.Consequently, we have As the sums in Eq. ( 5) are either equal to 0 or 1, resulting from Eq. ( 4), only two options are possible: If the edge is colored by c, the relation simplifies to otherwise, it equals to showing that (5) remains valid in both cases.
After applying the same approach to Eq. ( 2) and writing constraints through the entire area A, we obtain the following binary linear program: Despite the program (8) being formulated straightforwardly, it is hard to solve in this original form due to the (non-convex) integrality constraint (8e).Note that the problem is NP-complete because of its combinatorial nature [31].

Linear Programming Relaxation
In order to make the problem (8) easier to solve, let us relax the integrality constraint (8e) into a linear form Clearly, Eq. ( 9) allows for all configurations of (8e), but additionally also intermediate non-binary values.In addition, because the constraints in (8) are linear, the resulting approximation is convex and reads as

Semidefinite Programming Relaxation
The integrality constraint (8e) can be rewritten using the non-convex quadratic constraint which is satisfied only if x i,j,k = {0, 1}, being thus equivalent to (8e).This quadratic form does not simplify the solution, because the problem is still NP-complete, but we will use it for the derivation of a semidefinite programming relaxation.
Let us now substitute the binary variables x ∈ {0, 1} with y ∈ {−1, 1}, which requires Inserting this into Eq.(3) yields Consequently, we can write a non-convex quadratically-constrained optimization program in terms of y Let us further introduce a symmetric matrix Y ∈ S nt,y•nt,x•nt defined as The definition directly implies that any solution to the quadratically-constrained optimization problem ( 14 k∈T The only non-convex constraint (16f) can be relaxed, see, e.g., [32], into a convex form Y − yy T 0, (17) and based on the Schur complement lemma equivalently rewritten to Finally, the semidefinite programming relaxation of the binary linear program (8)

Extensions
Introduced formulations can include additional requirements for generated tilings.We list some of them below, but only in terms of x, as substitution of Eq. ( 12) into developed equations is obvious.

Tile-Based Boundary Conditions
There are four types of tile-based boundary conditions.First, we can enforce placement of the tile k ∈ T at position (i, j): Conversely, avoiding the tile k at (i, j) requires The requirement that the same tile is placed at both (i, j) and (p, q) can be written as Enforcing different tiles at (i, j) and (p, q) requires

Edge-Based Boundary Conditions
Starting from Eq. ( 5), we can also easily formulate edge-based boundary conditions.To constrain the north edge at (i, j) to the color c ∈ C, we write If the north edge at (i, j) has to differ from c, it holds that Equal color c of two edges, e.g., of the north edge at (i, j) and of the west edge at (p, q), is provided by Finally, different coloring of the north edge at (i, j) and the west edge at (p, q) is guaranteed through

Periodic Tiling
Periodic tiling is a tiling with equally colored both horizontal and also both vertical edges.Thus, periodic tiling secures tilability of the infinite plane.Building on the edge-based boundary conditions, a periodic tiling can be enforced by

. Tile Packing Problem
The tile packing problem, see, e.g., [30], consists in generation of a tiling with periodic boundary conditions, and each tile has to be placed exactly once.
In order to describe the problem in our framework, we just need to use Eq. ( 28) together with i∈H j∈W x i,j,k = 1, ∀k ∈ T . (29)

Objective Function
Obviously, the developed formulations can contain arbitrary convex objective function.For instance, we can minimize the maximum occurrences of tiles, compose the tiling such that it fits some pattern, etc.

The Branch and Bound Method
The whole integer design space, i.e., {0, 1} for the linear (10) and {−1, 1} for the semidefinite (19) approximation, can be described using a tree data structure, with each node corresponding to variables x, or y.However, as there is exactly one tile per position, it is advantageous to use the n t -ary tree representation instead of the binary one, so that the nodes correspond to the positions and branch into n t children, avoiding solution to infeasible programs placing multiple tiles at the same position.Solution to the developed approximations, (10) or (19), corresponds to exploring the root node of the tree.Unfortunately, because the approximations widen the feasible design space, their solution does not assure to solve the original problem (8) due to the design variables being allowed to take non-integer values.In order to overcome that, we branch some nodes of the tree, i.e., gradually fix all variables that correspond to the specific tile position to all possible combinations of integer values, n t in our case.
The relaxation is solved in each node, bounding the problem 1 , hence the name of the method branch and bound [33].The branching and bounding continues until the optimal (feasible) solution to (8) is found.If no such solution exists, the solution space is proven to be empty and the tiling generation problem infeasible.Without boundary conditions, infeasibility proves the tile set does not tile infinite plane.

Branching Rule
In our implementation, we firstly branch the most promising nodes, in which the integer infeasibility of the convex relaxation is minimal.For the linear relaxation, we have or for the semidefinite relaxation.If I inf = 0, a feasible solution to the original problem was found.

Variables Ordering
In each node to be branched, we have to define variables that will be fixed.It seems to be advantageous to fix the variables corresponding to the most difficult position, i.e., to the position with the highest integer infeasibility.In the case of the linear relaxation, the to-be-fixed position (i, j) requires max ∀i∈H,∀j∈W k∈T Similarly, for the semidefinite relaxation we write max ∀i∈H,∀j∈W

Examples
Building upon the previous sections, we implemented a custom branch and bound algorithm in C++.The linear relaxations are solved using Gurobi [34], the semidefinite programming relaxations are optimized by Mosek [35].As Gurobi contains a build-in state-ofthe-art branch and cut algorithm, usable only for the linear relaxations, we also measured its performance.Five tile sets were tested altogether: Jeandel's 11 tiles over 4 colors [16], Čulík's 13 tiles [15], Ammann's 16 tiles [12], Lauchli's 40 tiles [12], and Nurmi's set of 30 corner tiles [19].The tile sets were selected to sample different construction methods, tile set sizes, and numbers of colors used.Results of the benchmarks, ran on the Intel ® Core ™ i5-4210H processor, are summarized in Table 1.They indicate that the linear relaxation surpassed the semidefinite one in terms of speed and in the number of explored nodes.The latter results from the property of the simplex algorithm, used for the solution of linear relaxations, to terminate in vertices of the feasible space polytope.These are more likely integer-feasible than an interior point found by the interior-point algorithm, used for the solution of semidefinite programs.
Comparison of our branch and bound algorithm and Gurobi's build-in branch and cut, both using the linear relaxation, is even clearer.Our implementation lacks heuristics, parallelism, and generation of cuts, consequently being inefficient for tile sets with large number of degrees of freedom, such as the Čulík one.

Conclusions
In this contribution, we demonstrated that the tiling generation problem of finite-sized areas can be formulated as an integer optimization problem; and we also introduced its linear and semidefinite programming relaxations.Both the relaxations were employed in the branch and bound method and benchmarked on five tile sets.The results indicate that the linear relaxation is more suitable for practical applications.

8
) renders the matrix Y positive semidefinite, in our notation Y 0, and the rank of the matrix Y equal to 1.Moreover, the definition of y constrains all the elements in the main diagonal of Y equal to 1. Consequently, another equivalent formulation to the problem (,j,k + 1)[e k = c]

1
Bounding occurs only if an objective fuction is employed.

Table 1 .
Time demands and solved relaxations count in generation of aperiodic tilings.LP and SDP denote custom developed branch and bound method supplied with the linear and semidefinite relaxation, respectively; LP G stands for Gurobi cut-and-branch method and the linear relaxation.
East edge color of the tile k i Row iterator j Column iterator k, , m Tile iterators nc Number of colors in C n k North edge color of the tile k s k South edge color of the tile k n t,h Number of rows in H nt,w Number of columns in W x i,j,k Design variable, x i,j,k ∈ {0, 1} y i,j,k Design variable, y i,j,k ∈ {−1, 1} w k West edge color of the tile k I inf Integer infeasibility