Magic Squares

\[\begin{bmatrix} 🐯&😅&💫&🎩 \\ 🐰&🐇&🌙&🌟 \\ ✨&🎉&🎱&🃏 \\ 🔮&🔓&✅&🚀 \\ \end{bmatrix}\] The Magic Squares game is a simple yet engaging problem. There are many variations and the one I am interested in has the following rules

• 1. Use 4 by 4 matrix
• 2. Each cell can have a number from \(\{1,2,3,4\}\)
• 3. All rows should add up to 10
• 4. All columns should add up to 10
• 5. All diagonals should add up to 10

Here is an example of an invalid solution. While its rows add up to 10, its columns do not. So how one would go about finding a valid solution?

\(\begin{bmatrix} 1&2&3&4 \\ 1&2&3&4 \\ 1&2&3&4 \\ 1&2&3&4 \\ \end{bmatrix}\)

Brute Force Approach

Well, we can start to search by evaluating every possible 4x4 matrix until we find a valid solution. But how many matrices are there to evaluate? How do we generate them? Thankfully, a 4x4 matrix can be represented as a 16 digit number. Here is one, let’s call it \(s\), \(1111 \ 1111 \ 1111 \ 1111\). That would be a starting point. We can be certain that \(4444 \ 4444 \ 4444 \ 4444\) is the end and lets name it \(e\). Everything in between is a potential solution, so we have \(e-s\) which ends up being over 3 quadrillion candidates or 3 million billion.

I wrote functions to convert between 16 digit numbers and 4x4 matrices and check whether the solution is valid according to the game rules.

num_to_vec = function(n) {
  '
  Convert a number into a vector of its constituent digits.
  '
  vec = as.numeric(strsplit(as.character(n), "")[[1]])
  stopifnot(length(vec) == 16)
  return(vec)
}

vec_to_mat = function(v) {
  '
  Convert a 16 size vector to a 4x4 matrix.
  '
  if (typeof(v) == 'double'){v = num_to_vec(v)}
  stopifnot(length(v) == 16)
  return(matrix(v, ncol=4, byrow=TRUE))
}


check_solution = function(sol) {
  '
  Check if solution is valid.
  '
  # vector to matrix
  sol = vec_to_mat(sol)
  
  # diagonal positions
  pos_diag1 = c(1,6,11,16)
  pos_diag2 = c(4,7,10,13)
  
  # col and row sums
  d1 = sum(sol[pos_diag1]) == 10
  d2 = sum(sol[pos_diag2]) == 10
  
  col_sums = length(which(colSums(sol) != 10)) == 0
  row_sums = length(which(rowSums(sol) != 10)) == 0

  if(d1 & d2 & col_sums & row_sums) {
    print(sol)
    print('Solution is valid')
    return(TRUE)
  } 
  return(FALSE)
}

# define start and end 
s = 1111111111111111
e = 4444444444444444

# convert number to vector
num_to_vec(s)

[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

# convert vector to matrix and print
write_matex(vec_to_mat(s))

\(\begin{bmatrix} 1&1&1&1 \\ 1&1&1&1 \\ 1&1&1&1 \\ 1&1&1&1 \\ \end{bmatrix}\)

# check if matrix is a valid solution
check_solution(s)

[1] FALSE

Now, we are ready to start looping through the quadrillions of matrices. You will find out that this loop is unlikely to complete on a domestic computer. It will not even reach one valid solution before the computer runs out of memory.

# loop through all possible 4x4 matrices
candidates = s:e
for (c in candidates) {
  res = check_solution(c)
  if(res){print(c); break}
}

If we just try to load the last element in the candidates vector (which willl generate the whole vector), the R interpreter lets us know that we are short of just 12K GB of memory.

candidates[-1]

A generator with a yield statement would be my next step. Ultimately, this brute force approach is too naive, does not take advantage of any optimizations and it probably is better to implement in C or some other low level language. Do we have other ways to get just one valid solution with the help of R? Thankfully, we do.

Integer Linear Programming

ILP is a great tool to have in your toolkit. ILP enables us to find optimal solutions for a given a set of constraints and an objective function using linear optimization.

Below is a simple example. \(x_1\) and \(x_2\) are implicitly assumed to be integers. In some cases, they also may be set to binary, either a \(0\) or \(1\) .

\[\begin{align} \max \ x_1 + x_2 \\ \text{subject to} \notag -x_1 + x_2 \leq 12 \\ 3x_1 + 2x_2 = 5 \end{align}\]

The game rules I stated in the beginning define linear relationships between the elements of the 4x4 matrix. We can translate those rules into the constraints and the magic squares game into an Integer Linear Program. So let’s do that. Let’s start with the objective function.

We already know that we have 16 elements in a 4x4 matrix, for each element we may use one of the four numbers \(\{1,2,3,4\}\). This means we have \(16*4 = 64\) binary variables. Our objective function will be the sum \(x_1,x_2,...,x_n\) where \(n\) is \(64\). Great. Normally we would say we want to maximize the objective function. Here we know that we need exactly 16 elements for any valid solution out of 64. So we can just set the objective function to be equal to 16. I will write down our ILP in notation and explain each constraint.

\[\begin{align} \sum_{n=1}^{64} x_n = 16 \\ \text{subject to} \notag \sum diag(x_n) = 10 \ & (1) \\ \sum col(x_n) = 10 \ & (2) \\ \sum row(x_n) = 10 \ &(3) \\ \sum cell(x_n) = 1 \ & (4) \\ \sum digit(x_n) = 4 \ & (5) \end{align}\]

I am certainly taking liberty in the way I use notation here. (1), (2), (3) are straightforward and make sure that \(x\)’s that correspond to rows, columns and diagonals sum up to 10. (4) makes sure that \(x\)’s that correspond to the same cell in the matrix sum up to 1. Indeed, each cell in the matrix can only have one number from four available and thats why only one corresponding \(x\) can be non zero. (5) ensures that all four available numbers are used at least once. Implicitly, the \(x\)’s are binary.

So we actually ended up having a system of equalities. Let’s solve this system using the lpsolve package.

lpsolve will need the objective function, the constraint matrix, the equality/inequality direction, and finally the right hand side of the equalities/inequalities. I code everything below.

# Helper functions to build constraints and process solutions

position = function(el) {
  '
  Get vector of positions of elements in a ncols(el)*4 
  matrix that correspond to positions in the original 4*4 matrix.
  '
  x1 = el*4-4+1
  x2 = el*4
 
  return(x1:x2)
}

constraint = function(el_pos,val=1, raw_pos = FALSE) {
  '
  Create constraint vector with 64 elements given element positions 
  '
  base_vec = rep(0,4*16)
  if(!raw_pos) {poss = c(sapply(el_pos, position))}
  else {poss = el_pos}
  
  base_vec[poss] = rep(val,length(poss)/length(val))
  return(base_vec)
  
}

# define coefficients in objective function
obj = rep(1,4*16)

# set up constraints base vector
cons_base = rep(0,4*16)

# diagonal sums are 10
pos_diag1 = c(1,6,11,16)
pos_diag2 = c(4,7,10,13)

con_diag1 = constraint(pos_diag1,val=1:4)
con_diag2 = constraint(pos_diag2,1:4)

diag_cons = rbind(con_diag1,con_diag2)
rownames(diag_cons) = rep('diag_sums',2)

# col sums are 10
col_cons = c()
col_positions = matrix(seq(1,16),nrow = 4)
for(row in 1:nrow(col_positions)) {

  col_cons = rbind(col_cons, constraint(col_positions[row,],1:4))
}
rownames(col_cons) = rep('col_sums',4)

# row sums are 10
row_cons = c()
row_positions = matrix(seq(1,16), ncol = 4,byrow=TRUE)
for(row in 1:nrow(row_positions)) {
  row_cons = rbind(row_cons,constraint(row_positions[row,],1:4))
}
rownames(row_cons) = rep('row_sums',4)

# each cell can have only one value out of 4 available
only_one_val = t(sapply(1:16,constraint))

# each number must appear 4 times
appear4 = rbind(constraint(seq(1,64,by=4), raw_pos=TRUE),
                constraint(seq(2,64,by=4), raw_pos=TRUE),
                constraint(seq(3,64,by=4), raw_pos=TRUE),
                constraint(seq(4,64,by=4), raw_pos=TRUE))
  

# create constraint matrix
cons = rbind(diag_cons,col_cons,row_cons, only_one_val, appear4)

The cons matrix has 64 columns and a constraint for each game rule we outlined. Let’s define relationship direction and right hand sides.

# Set inequality/equality signs
f.dir <- c(rep("=",nrow(diag_cons)),
           rep("=",nrow(col_cons)),
           rep("=",nrow(row_cons)),
           rep("=",nrow(only_one_val)),
           rep("=",nrow(appear4))
           )

# Set right hand side coefficients
f.rhs <- c(rep(10,nrow(diag_cons)),
           rep(10,nrow(col_cons)),
           rep(10,nrow(row_cons)),
           rep(1,nrow(only_one_val)),
           rep(4,nrow(appear4)))

Finally, we are ready to solve this ILP and see if can obtain at least one valid solution without runnning out of memory or searching through a quadrillion numbers.

library(lpSolve)

# Solve ILP
model = lp("max", obj, cons, f.dir, f.rhs, all.bin=T,use.rw=TRUE)
print(model)

Success: the objective function is 16 

model$solution

 [1] 0 0 0 1 1 0 0 0 0 0 0 1 1 0 0 0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 0 1 0 0 0 0 0 1 1 0 0 0 0 0 1
[48] 0 0 0 1 0 0 1 0 0 0 0 1 0 0 1 0 0

Ok, great. Seems we got a solution quite quickly, almost instantly. Looks confusing. Above vector shows us which 16 of 64 \(x\) variables are activated and can be included in the solution matrix. I define a function to translate the solution vector into a matrix. Let’s check!

collect_solutions = function(sol) {
  '
  Print out the solution matrices.
  '
  sol = sol[1:(as.integer(length(sol)/64)*64)]
  base_vec = rep(1:4,16)
  sol = matrix(sol,ncol=64, byrow = TRUE)
  solutions = c()
  for(s in 1:nrow(sol)){
    sol_mat = matrix(base_vec[which(sol[s,]>0)], ncol = 4)
    solutions = rbind(solutions,c(sol_mat))
  }
  
  return(solutions)
  
}

sol = collect_solutions(model$solution)
sol # the vector has 16 elements. Great!

     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
[1,]    4    1    4    1    1    3    2    4    2     4     1     3     3     2     3     2

check_solution(sol)

     [,1] [,2] [,3] [,4]
[1,]    4    1    4    1
[2,]    1    3    2    4
[3,]    2    4    1    3
[4,]    3    2    3    2
[1] "Solution is valid"
[1] TRUE

After checking our solution according to the game rules we see it is valid. We basically created a solver for a variation of the magic squares game using ILP and lpsolve. In general, any problem with linear constraints and an objective can be solved using Linear Optimization. Hope this illustration was informative to you.

\[\begin{bmatrix} 🕛&🕧&🕐&🕜 \\ 🕑&🕝&🕒&🕞 \\ 🕓&🕟&🕔&🕠 \\ 🕕&🕡&🕖&🍀 \\ \end{bmatrix}\]