Fitting a pairwise binary Markov random field by regularized moment matching and sampling the resulting Boltzmann distribution with parallel-tempered Markov chain Monte Carlo(to get better at Pokémon)
I'll walk you through understanding how the physics of magnets can be used to build competitive Pokemon teams (or understand voting patterns of the United States Supreme Court). The math might be a little intimidating, but this page is filled with interactive widgets to help you understand exactly what each equation is doing and how everything fits together.
A any sort of magnet can be thought of as composed of individual tiny magnets, each with it's own magnetic field, called a "spin". These spins can point either up or down, also called +1 and -1. When the spins all align, their magnetic fields combine to create the strong magnetic effect we can see macroscopically. The more aligned the spins are, the lower the energy of the system, and the stronger the magnet. Since nature prefers low energy states, these spins will align by themselves and create a macroscopic magnet.
However, spins can be randomly flipped by thermal fluctuations. Flips into lower energy states are more likely, but if you add enough heat then the spins will flip so often that they can't maintain their alignment. This is when the magnet "demagnetizes". The temperature this happens at is called the Curie point.
The Ising model is a way to describe that same behavior mathematically. We create a grid of cells, called a lattice, with each cell containing a "spin" that can point either up or down.
The energy of a configuration of spins with value is given by the formula
and the probability of seeing configuration is . This is called the Boltzmann distribution (remember that name, it'll come up again).
In basic terms, it means that the energy of the system is lower when more neighbors are aligned, and that lower energy states are more likely to occur, especially at low temperatures. This 2d lattice Ising model has been well-studied, and many aspects of its behavior have been solved analytically (without needing to simulate it) (Onsager, 1944). For example, the dotted red line in the plot below shows the expected magnetization for any temperature below the critical temperature (~2.269), where the system transitions between ordered and disordered phases.
The widget below lets you run a simulation of the Ising model and see how it behaves at high and low temperatures. This temperature slider does the exact same thing as the one in the Team Completer, as I'll explain further down. Take a moment to play with the simulation. See if you can use it to understand why magnets will sometimes flip their poles when they heat up and cool back down.
The 2d lattice model works well for explaining magnets, but we can also think of it as a special case of a more general model. The spins don't need to be in a grid, so we can represent any topology we want as a network of nodes (spins) and edges (couplings). We can also say that each spin has a bias , representing a preference for how it wants to spin (for example, caused by an external magnetic field).
The energy formula becomes:
If you plug in 0 for and 1 for , you'll get the exact energy formula we used for the lattice version.
Couplings can also be negative, meaning that those spins prefer to be opposite, rather than aligned. Also, from now on we'll think of spins as on or off (1 or 0) instead of +1 or -1. It makes more sense for what's coming up, and it doesn't change the math.
The figure below is a randomly generated network representing a more general Ising model. While math can predict the Boltzmann distribution of the simple 2d lattice analytically, that isn't yet possible for arbitrary graphs. Try a few random graphs at different temperatures to see how even a tiny system can behave completely unpredictably.
I've let you play with a few Ising simulations now, but I haven't actually explained how to go from mathematical equations to simulation. The key is the Markov Chain Monte Carlo (MCMC) method. The most basic version of MCMC, called the Metropolis algorithm (Metropolis et al., 1953), just looks at one spin at a time and decides whether to flip it based on one rule:
In other words, if flipping the spin would lower the energy of the system, then we flip it. If not, we still might flip it with a probability that increases as the energy difference gets smaller or as the temperature gets higher. Given infinite time, this process will find itself in configuration at the exact rate predicted by the Boltzmann distribution.
Click "Propose a flip" to set up a single Metropolis step.
Before we go further, it's important to understand the concept of an energy landscape. If you're already comfortable with the concept of energy landscapes then skip right on ahead to the next section. If not, I'll give you a quick primer here so that you can understand some of the more mind-bending concepts coming up.
An energy landscape is a way of visualizing the energy function that defines a system. Each point in the landscape corresponds to a particular configuration of the system, and the height of the landscape at that point corresponds to the energy of that configuration. A "configuration" here means a specific arrangement of all the spins, so the small 32x32 lattice you played with earlier has possible configurations.
Imagine a large mountain range: Each of those configurations is a different spot in that mountain range, and the altitude at each point is the energy of that configuration. Valleys represent low-energy configurations (like a fully aligned magnet) and peaks represent high-energy configurations. Systems, including the Ising model, tend to "roll downhill" in this energy landscape as they look for low-energy configurations. The Metropolis algorithm can be thought of as a hiker exploring this landscape: every step, he points in a random direction. If it's downhill, he always takes it. If it's uphill, he might still take it, but he's less likely to if it's a steep climb.
(Technically, the landscape in our case would be in a 1024-dimensional space, but picturing a 2d landscape (with the 3rd dimension representing energy) is an easy mental image and the math is exactly the same.)
Energy landscapes show up all over the place in science. The "mountain range" analogy is so common that it's common to use geographical terminology to describe what's going on. When I use terms like "valley", "ridge", or "basin" in the next few sections, I'm referring to features of the energy landscape. There are some simulations below that will make this all make a lot more sense.
If you got upset when I said earlier that Metropolis will sample the Boltzmann distribution "given infinite time", good instincts. The random walk through the energy landscape can get stuck in a specific basin for longer than the duration of a finite simulation. You can see this in action by running the simulation below.
Parallel tempering (Hukushima & Nemoto, 1996) runs chains at temperatures . The hottest chain roams over the whole landscape; the coldest stays concentrated in whichever well it currently occupies. Periodically, the algorithm proposes that adjacent chains swap configurations, with an acceptance probability that preserves the Boltzmann distribution:
The hotter chains explore the full landscape, and periodically swap with colder chains for detailed exploration. This gives a much more thorough sampling of the entire landscape relative to the single chain Metropolis version. Check for yourself by switching back and forth between single-chain and PT modes in the simulation.
| steps | 0 |
|---|---|
| cold chain in left well | 0.0% |
| true (by symmetry) | 50.0% |
Parallel tempering gives reliable results but it's computationally expensive, requiring thousands of sweeps across multiple chains. Mean field theory gives us an alternative which is much simpler to compute, at the cost of being only an approximation. Rather than attempting to sample the full Boltzmann distribution, MF directly estimates the mean spin values by iteratively applying the MF update equation:
On the energy landscape, this is gradient descent: just roll directly downhill. The marker below converges to whichever local minimum it starts nearest and stops. No swap moves, no temperature, no escape. This is similar to (though not exactly the same as) running Metropolis or Parallel Tempering at .
We also use it as the starting point for the Team Completer if you don't use full sampling. We were talking about Pokemon, remember? We'll get back to Pokemon and how this all ties together soon, but first I need to take one more detour and talk about the United States Supreme Court from 1994-2005.
| iterations | 0 |
|---|---|
| |gradient| | 1.42e+0 |
| status | descending |
| settled in | — well |
Loading SCOTUS data…
Alright, it's finally time to actually talk about Pokémon! By now you should have a good understanding of the methods, so all that's left is to connect the dots. Starting with the data: thanks to the amazing people over at Limitless VGC, we have a huge dataset of real competitive teams that players have brought to tournaments. Limiting to the current regulation, doubles tournaments, and only events with at least 64 participants, we get about 4,000 teams. These are the observations the model is fit to.
In the SCOTUS example there were 9 spins to fit, one per justice. The current VGC regulation, Pokémon Champions Regulation M-C, has roughly 200 species and formes and about 100 held items. Modeling species alone is around 200 spins; every Species @ Item combination is closer to 20,000. With one parameter per spin plus one per pair of spins (for the coupling), that full model needs to fit about 200 million parameters (nearly a gigabyte), all of which still has to run in the browser! So we keep only the spins that appear in 5 or more teams, bringing the model down to 318 spins and about 50,000 parameters, about 1% of the full pairwise space. Much nicer.
The largest difference between the Pokémon and SCOTUS/Ising cases is that in those earlier examples every combination of spins was valid. Any spin could be up or down regardless of the others, and any justice could vote liberal or conservative regardless of the others. In Pokémon that's not true: only six Pokémon can be on a team at once, and duplicate species or items are forbidden. These states aren't just unlikely, they're literally impossible. To enforce this, we make a small modification to the sampler that still produces the correct Boltzmann distribution over the valid states. Instead of considering one spin at a time and proposing to flip it, we start with six spins on and consider two spins at a time, one on and one off. The proposal is to swap their states, keeping exactly six on at all times. If a swap would violate the species or item constraints it is rejected outright (without this rule such states would rarely appear anyway since the model assigns them very high energy, but it's good to enforce it explicitly).
Because a feature here is a pair, a species together with its held item, the real sampler splits that swap into two kinds of move. One proposes exchanging a species on the team for one off it, scoring the proposal against every item that species could plausibly hold rather than a single build, then drawing the item it actually arrives with. The other leaves the six species alone and just rerolls one of their items. Picking randomly between the two on each step explores which Pokémon belong on the team and which items they should hold as separate questions, which mixes far better than treating every species and item combination as its own indivisible spin.
The figure below shows how the couplings in the fitted model connect the top species together. Use the slider to hide the weaker couplings and watch the structure of the strongest relationships emerge.
And that's it! Hopefully this page has given you a feel for how k2dex builds and analyzes teams using statistical physics. Building this system and this page has been a ton of fun for me, and I hope you find it both interesting as an application of statistical physics and useful as a teambuilding resource!
This isn't a research paper, but the methods here are well-established and worth crediting. These are the works that the model and this explainer build on, in roughly the order their ideas appear above.