MIP: Wood puzzle

I found a wooden puzzle in my mother's bookshelf. This is all there is to it:

All pieces can be moved freely. A single piece:

Each piece has one corner each in the colors  Red, Green, Blue, and Yellow. The goal is to match colors of adjacent corners. Solving this becomes a neat programming exercise.

Brute Force

How long would it take to just enumerate all solutions? There are 8 pieces and 8 locations, which gives us 8! = 40320 permutations. However, the puzzle can be rotated 180 degrees and be identical, so it's just 20160 permutations of the pieces. The individual pieces can also be turned in 4 different directions independently of each other. So for each permutation, we have 4^8 = 2^16 = 65536 turns. In total we have 8! / 2 * 2^16 ~ 1.3 billion solutions. So it takes about the same amount of time for the computer to solve it by brute force, as for a person to solve it manually.

MIP

The problem is naturally modeled as a constraint problem. For MIP, one can use one binary variable for each combination of location, piece, and turn. In total, 8*8*4 = 256 variables. There are 22 "conflict edges" where colors from two different pieces must match in an adjacent corner.

The most numerous constraints, however, are the structural constraints. We must ensure that each piece gets exactly one location, and vice versa. Each piece must also get exactly one turn. These constraints are all on the same form: a sum of a set of binary variables must equal 1. These are also called symmetric constraints.

For each corner of a location, we must know its color, as a function of the decision variables. This is necessary so that we can express the color matching constraints in terms of the decision variables. For each corner of a location (dot) and each color, I create one auxiliary variable. In total, that becomes 8*4*4 = 128 dot-by-color variables.

Pattern Picking constraints

The dot-by-color variables are bound to the decision variables using pattern picking constraints. For each decision variable X and each dot-by-color variable Y, make the constraint:

If X dictates the same color for Y's dot as Y does, then:

X <= Y

Otherwise:

X + Y <= 1

How does this work? If X is 0, then Y can be either 0 or 1 without violating the constraints. However, if X is 1, then Y is forced to be 1 or 0 respectively, depending on which constraint is active. Note that we use the same Y variables for all decision variables, but they are cleverly bound together to all of them without risk for infeasible solutions.

Extra constraints for efficiency

It is implicit from the constraints on the pieces that for a given location, there will be exactly one dot per color, and one color per dot. This is something we can tell the solver beforehand, as constraints. This makes the solution go almost 4 times faster.

Result

The solution takes about 8 seconds for find with CBC, so it is not trivial for it.

Conclusion

Most constraints are for making the model internally self-consistent. These types of constraints arise from things that humans take for granted, such as the pieces being non-reusable. These internal constraints must be uncovered during the modeling process, and take a lot of time.

The pattern picking constraints are quite neat and abstract, and could possibly be packaged as part of a solver interface.

We can do some of the work for the solver by adding constraints that we can logically infer from the problem. This can be important for performance optimizing MIP models.

Playing with Word2Vec

In a lecture by Yann LeCunn, he mentions that word2vec embeddings can be used to learn analogy properties. I want to see if any other natural language concepts are modeled by word2vec's.  

A word2vec embedding is a vector on the order of 100 elements that represents a word. The embedding is the hidden layer of a network that for each pair of words predicts the probability that they appear close to each other in a large corpus. An example of an analogy:

London is to England what Paris is to France.

When people have looked at word2vec's, they have found that the following property holds approximately true:

w2v(london) - w2v(england) = w2v(paris) - w2v(france)

This difference is also the same for most capital-country pairs. That's quite interesting! Let's see if we can find some other interesting properties.

I use a dataset of 43000 words from Kaggle. Many of the words have leading upper case, while not being proper names. After converting everything to lower case and removing duplicates (keeping only the first to appear in the list), there remains 38000 words. 

Relation
First, wouldn't it be great if this was true?

w2v(london) - w2v(england) = w2v(capital)

Looking at the 10 closest elements of london - england:

paris: 0.345
tokyo: 0.321
brussels: 0.315
mayfair: 0.314
uptown: 0.303
nairobi: 0.302
jakarta: 0.302
budapest: 0.299
amsterdam: 0.294

OK, that didn't work. It doesn't work to take england - london either. The similarity to "capital" is only 0.15.

Hypernym
A hypernym of a set of words is a word that describes the set as a whole. For example:

Cutlery is a hypernym of knife, fork, and spoon

Color is a hypernym of red, green, and blue

Can we shake out the hypernym with word vectors?

hypernym(words) = mean([w2v(word) for word in words])

hypernym(["red", "green", "blue"])
yellow: 0.459
purple: 0.440
orange: 0.436
pink: 0.433
brown: 0.418
purple: 0.418
white: 0.407
gray: 0.396
colored: 0.387
crimson: 0.380
maroon: 0.372
pink: 0.370
color: 0.355

Not too bad! Color shows up among the top answers.

w2v.hypernym(["mercedes", "jeep", "ford"])  

1 mercedes: 0.615
2 ford: 0.535
3 jeep: 0.514
4 car: 0.511
5 jeep: 0.505
6 sedan: 0.490
7 jaguar: 0.470
8 buick: 0.468
9 vehicle: 0.458
10 mustang: 0.448

Great! Car is the top answer except for the three included words. 

w2v.hypernym(["bed", "table", "chair"])                                                                                                                                           

1 bed: 0.499

2 chair: 0.474

3 table: 0.410

4 sofa: 0.384

5 couch: 0.378

6 divan: 0.359

7 beds: 0.349

8 chairs: 0.334

9 footstool: 0.333

10 daybed: 0.328

I wanted to see 'furniture' here, but no dice.

Reduce: 'country'
What if we subtract more than one vector from another? Can we peel back layers of meaning this way? Example:

The word 'country' has several meanings. It can be a music genre, it can be a synonym for 'nation', and it can refer to land outside of cities. Remove the 'land' and 'nation' context, we should be thinking of music.

For the implementation of this, I tried two approaches: simply subtracting away vectors, and also projecting away components. They both worked about as well. The word vector is normalized after each reduction, so the scale is relevant throughout.

w2v.reduce("country", [])                                                                                                                                                         

2 nation: 0.724

3 world: 0.598

4 globe: 0.514

5 america: 0.486

6 countries: 0.482

7 national: 0.478

8 abroad: 0.472

9 republic: 0.453

10 europe: 0.453

Without context, the strongest associations to country is in the sense of 'nation'.

w2v.reduce("country", ["nation"])                                                                                                                                                 

2 abroad: 0.445

3 countryside: 0.374

4 homeland: 0.363

5 border: 0.347

6 province: 0.347

7 villages: 0.346

8 rumanians: 0.336

9 europe: 0.335

10 thailand: 0.335

And with a heap of reductions we can indeed force it to think only of music!

w2v.reduce("country", ["land", "nation", "rural", "abroad", "province", "europe"])                                                                                                

2 genres: 0.211

3 carreer: 0.207

4 singer: 0.206

5 musician: 0.177

6 catatonia: 0.175

7 promoter: 0.173

8 roped: 0.169

9 duet: 0.166

10 crooning: 0.165

11 vocalist: 0.164

Reduce: 'house'

Can we do the same for 'house'? Also a music genre, with more than one meaning.

w2v.reduce("house", [])                                                                                                                                                           

2 senate: 0.702

3 bill: 0.541

4 appropriations: 0.520

5 congress: 0.511

6 commons: 0.487

7 congressional: 0.486

8 senators: 0.483

9 lawmakers: 0.473

10 conferees: 0.447

I make a rule to always reduce the top word.

w2v.reduce("house", ["senate"])                                                                                                                                                   

2 manor: 0.461

3 sanctuary: 0.414

4 lodge: 0.409

5 lounge: 0.400

6 gardens: 0.388

7 club: 0.388

8 palace: 0.381

9 gate: 0.377

10 factory: 0.375

... repeating this a number of times ...

w2v.reduce("house", ["senate", "manor", "factory", "lords", "society", "rowdy", "sanctuary"])                                                                                     

2 ways: 0.155

3 true: 0.133

4 diet: 0.132

5 networks: 0.122

6 inset: 0.120

7 r: 0.116

8 feat: 0.114

9 signature: 0.113

10 crestfallen: 0.110

Yes! At least we get 'feat' (short for featuring) and 'true' in there.

Alignment
Can the embeddings be used to sort things? Example:

Since the sun is bigger than a car, 'sun' should be closer to 'big' than 'car' is. 

This is admittedly quite a long shot for such a simple model, and it turns out to be very wrong:

w2v.alignment("big", ["sun", "planet", "car", "dog", "ant"])                                                                                   

1 car: 0.120

2 ant: 0.100

3 dog: 0.099

4 planet: 0.079

5 sun: 0.022

w2v.alignment("healthy", ["carrot", "burger", "wine", "juice", "cake"])                                                                                         

1 juice: 0.221
2 carrot: 0.161
3 cake: 0.084
4 burger: 0.054
5 wine: 0.025

Probably what we are seeing here is mostly how often people talk about these objects in the context of being "big" or "healthy".

Conclusion

Word embeddings can yield surprisingly relevant results with very simple models. However, it is important to know that they have been generated using a simple optimization that only looks at local context, and has no real way of modelling that context. This becomes apparent when trying to make it perform more abstract tasks.

Searching expressions

During a coffee break at work, the following brainteaser came up:

0 0 0 = 6
1 1 1 = 6
2 2 2 = 6
3 3 3 = 6
4 4 4 = 6
5 5 5 = 6
6 6 6 = 6
7 7 7 = 6
8 8 8 = 6
9 9 9 = 6
10 10 10 = 6

For each row, one is supposed to fill in the operators that make the expression true. One is not allowed to add more digits, of course. The allowed operators are the arithmetic operators, exponentiation, square root, and factorial. For example, we have:

(1 + 1 + 1)! = 6
3 * 3 - 3 = 6

We solved 0 through 10, but struggled with 11. We also wondered whether there were some interesting solutions for 0-10 that we had missed.

I have been intrigued by symbolic programming and automated theorem proving since I first read about it. This seemed like a simple enough problem to start exploring those subjects. The first obstacle, which proved to be the most difficult one in the end, was conceptual blocks on my part. From proving mathematical theorems in university, I'm used to feeling that at each step in the derivation, there is an "infinite" number of possible moves, if only one is creative enough.

Tokens. x, y (numbers)
Unary operator. U: x -> y
Merge operator. M: x, y -> z

Example
State: [x, y, z]
Applied operation: M1 ('add'), index 0
Results: [M0(x, y), z] = [x+y, z]

Searching through expressions
However, if we want to do a systematic search in the computer, we need to be able to enumerate the possible moves. In our case, we have restricted ourselves to a finite number of operators, and there is a finite number of tokens to apply them to. So the number of moves in each step is not only enumerable, but finite. My idea for looking for an answer here is to simply do exhaustive search, using BFS (breadth-first search). Breadth-first search should help us find one valid solution quickly (if we assume that there is at least one relatively simple solution). However, if we want to find all solutions, it doesn't really matter which search algorithm is used, as we will have to go through all states anyway.

The really important thing is termination conditions. I pruned all states which had at least one token which was not a real integer. I also pruned all states with an integer larger than 1000.

Implementation
The search can be done very conveniently using networkx. The states are represented as nodes (the key is a tuple of the tokens). The edges are directed, and each edge has a feature with a representation of the operation. This help verify correctness, and display the result.

Another question is how to return solutions. First, I tried to return all simple paths from the starting stated to the goal state. For x=9, this gave over 17000 solutions. Clearly a lot of them were the same up to permutation. I decided to return one solution per immediate predecessor to the target state. For each immediate predecessor, I return the shortest path in the search graph. This brought down the number of solutions to no more than a few 10's at most.

Results
The resulting successful computations can be plotted using graphviz via pydot. Installed through:

apt-get install graphviz
pip3 install pydot




Some examples. The _index suffixes are there because pydot doesn't accept custom labels.

The only solution found for 1.

One of two solutions found for 8.

One of 22 solutions found for 9.

The searcher found many solutions that I hadn't thought about, however most were not very interesting. I ran it for numbers up to 1000, and sadly there was not a lot of surprises. Going back to the problem we had in the beginning, it did not find a solution for 11. The workable numbers up to 1000 are typically very even numbers. The highest workable I found was 729, which goes through 81:

Path from 729. which goes through [729, 729, 729]->[27, 27, 27]->81->9->3->6.

The highest workable prime found is 37, which goes through 36:

Conclusion
This was a very fun exercise and it yielded successful results quickly, even if it was not so surprising. It should be noted that we're relying on numerical answers here as opposed to having an exact representation of numbers. This works well for integers, but wouldn't generalize to reals.

Convex Hull: correlated points

Previous:

• Convex Hull

In the previous post, I investigated the convex hull of random point clouds. I found that as we add more dimensions while keeping the number of points the same, the expected volume decreases sharply. One way to think of this is that in high dimensions, almost all points in a random point cloud will lie on the edge of the convex hull, very few will be "surrounded" by the other points. 

In this post, I want to see what happens if the points are not random, but correlated to each other. Why might we think so? Well, strongly correlated systems have in a sense fewer degrees of freedom, so we expect them to behave like lower dimensional problems. Like this:

Correlated points in R^2. Out of the 100, 90% of the points are in the interior of the convex hull of the other points. These are marked on blue. The points on the edge of the convex hull are marked in red.

How to sample correlated points in R^N

Generate a square matrix A with uniform random elements, each with expectation 0. For each point in the point cloud, generate N gaussian random numbers, and apply A to these. The covariance matrix will be AA^T. 

2 dimensions

Making 100 random trials gives that the average k-fold inclusion is 91%, compared to 88% for just random point clouds. So the difference is not big, not even statistically significant. The empiric variance of an estimate probability p in a binary distribution is n*p*(1-p), where n is the number of trials. With n=100 and p=0.9 in our case, this comes out to about 8. So the standard deviation is about 2.8, pretty much the same as the difference between the sets. We might want to apply a real test here, something that respects the binary distribution. But that's silly; if we wanted something more accurate we should just make more tests. 

10 dimensions

Doing this in 10 dimensions give an estimate of 5.6% for the independent point clouds, and 6.0% for the correlated point clouds. Once again, not significant. 

Conclusion

When we observe high-dimensional systems, almost all observations will be "extreme" in some direction, even if the system is structured. 

Sampling: simple stable Autoregressive

I want to try out the convex hull inclusivity tester on data that is less random. One way to get such data is to sample a stable ARMA process. In fourier space, an ARMA can be expressed as a quotient of two polynomials:



P contains the poles, and Q contains the zeroes, or nills. In this post, I use only Q = 1.
The criterion for stability is that all roots of P is in the unit circle. To get a uniform sampling from the unit circle, we can sample an angle from a uniform distribution, and an absolute value from a triangle distribution (since the density must increase linearly with the absolute value). The code:

This output format fits perfectly into the ArmaProcess object from statsmodels.

Random stable AR process from 5 pairs of poles.

Convex Hull

I want to find the convex hull of some set of points in  (a euclidean space of arbitrary dimension). I am hoping that this will be an interesting programming exercise, and that it will provide some insight into the properties of higher dimensions.

Definition

The Convex Hull of a set S is the smallest convex set that contains all elements of S. So much for the mathematical definition. Intuitively, we can think of the convex hull as the shape of a rubber band that encloses all points in S. 

Example of the convex hull of a set of points. 

A Simple Test

When starting out, I thought about explicitly constructing the convex hull. This is indeed something we will get to later. By explicitly constructing the hull, I mean finding the linear subspaces that constitute it. In the example above, that would imply finding the line equations of the red lines. However, it struck me that this is not necessary if all we want is to test whether some given point is inside the convex hull or not. There is a way to do it by solving a single linear programming problem. We use the following fact:

A vector u is in the convex hull of a set of points V if u can be expressed as a convex linear combination of the vectors in V. That is:

such that:

.

This can be readily formulated as a linear program, where u and V are given, and the lambdas are variables. The code:

This works beautifully in 2D:

Convex hull for the set of points (0, 0), (1, 0), (0, 1), marked with red dots. Sampled in the integers. 

Convex hull (red) of random points (green) sampled in the integers. 

Volume of Convex Hull using Monte Carlo

One question I had was: "suppose we take 10 random points sampled from the unit cube, what is the expected volume of their convex hull?". This requires a nested simulation. On the highest level, we need to sample a lot of 10-point sets and measure their volumes. And for each volume that we want to measure, we need to sample a lot of points to get an accurate Monte Carlo estimate of the volume. Some initial tests:

dim | volume: 10 points | volume: 100 points
2D  | 0.42              | 0.88
3D  | 0.14              | 0.69
4D  | 0.035             | 0.45

What can we say from this? Not much more than that the average coverage decreases as the number of dimensions increase. 

Monitoring Neural Network convergence with Multidimensional Scaling

I want to better understand what happens when a neural network is trained. One way is to analyze how the weights change. Since there is typically over 100 weights, we need to look at a low dimensional extract if we want a comprehensive view. I want to trace the search through weight space like a path, so we need to reduce to 2 or 3 dimensions.

An interesting approach to this extreme dimensionality reduction is Multidimensional Scaling (MDS). MDS uses a matrix of pairwise distances between the elements x_i to find low dimensional real valued z_i such that


is minimized, where d_ij is the distance between elements x_i and x_j, for some given metric. I used euclidean distance for the metric.

The problem

I apply this to a simple classification: given a 2-dimensional input, classify whether the point is in the unit circle. I use a neural network with 1 hidden layer of 30 nodes. Therefore, there are (2+1)*30 weights in the first layer, and (30+1)*1 weights in the output layer. (The +1 terms are the biases). This comes out to 121 weights. So we are reducing from 121 dimensions to 2, while trying to optimally preserve pairwise distances.

The search paths

Unsurprising: the search takes smaller and smaller steps as the training reaches a better state.
Surprising: the optimization appears to be just a straight march towards the goal. 

Looking at just the end of the search path. As the training converges, the search takes smaller and more erratic steps. But it does not walk in circles (yet, at this point). 

How can this be improved?

- Can look at each layer by itself

- Show MDS fit score, perhaps color code / size code the plot by uncertainty