Code for the Dark Worlds competition

Here’s the code for my entry in the Dark Worlds competition.

- Run ‘test_on_train.m’ to test the algorithm on the training data.
The settings in this file are such that you will get a good result in about 15 min on a quad-core PC.

- Run ‘predict_leaderboard.m’ to reproduce my winning submission for the competition. This requires a lot more running time for a relatively modest gain in accuracy.

Get the code.

Observing Dark Worlds

Kaggle recently ran another great competition, which I was very fortunate to win. The goal of this competition: detect clouds of dark matter floating around the universe through their effect on the light emitted by background galaxies.

From the competition website:

There is more to the Universe than meets the eye. Out in the cosmos exists a form of matter that outnumbers the stuff we can see by almost 7 to 1, and we don’t know what it is. What we do know is that it does not emit or absorb light, so we call it Dark Matter. Such a vast amount of aggregated matter does not go unnoticed. In fact we observe that this stuff aggregates and forms massive structures called Dark Matter Halos. Although dark, it warps and bends spacetime such that any light from a background galaxy which passes close to the Dark Matter will have its path altered and changed. This bending causes the galaxy to appear as an ellipse in the sky.

The task is then to use this “bending of light” to estimate where in the sky this dark matter is located.

Dark Matter bending the light from a background galaxy.

Dark Matter bending the light from a background galaxy.

Although the description makes this sound like a physics problem, it is really one of statistics: given the noisy data (the elliptical galaxies) recover the model and parameters (position and mass of the dark matter) that generated them. After recovering these dark matter halos, their positions could then be uploaded to the Kaggle website where a complicated loss function was used to calculate the accuracy of our estimates.

Bayesian analysis provided the winning recipe for solving this problem:

  1. Construct a prior distribution for the halo positions p(x), i.e. formulate our expectations about the halo positions before looking at the data.
  2. Construct a probabilistic model for the data (observed ellipticities of the galaxies) given the positions of the dark matter halos: p(e|x).
  3. Use Bayes’ rule to get the posterior distribution of the halo positions: p(x|e) = p(e|x)p(x)/p(e), i.e. use to the data to guess where the dark matter halos might be.
  4. Minimize the expected loss with respect to the posterior distribution over the predictions for the halo positions: \hat{x} = \arg\min_{\text{prediction}} \mathbb{E}_{p(x|e)} L(\text{prediction},x), i.e. tune our predictions to be as good as possible for the given error metric.

For step 1. I simply assumed that the dark matter halos were distributed uniformly at random across the sky. Step 2 is more complicated. Fortunately the competition organizers provided us with a set of training skies for which the positions of the dark matter halos was known, as well as a summary of the physics behind it all. After reading through the tutorials and forum posts it became clear that the following model should be reasonable:

p(e_{i}|x) = N(\sum_{j = \text{all halos}} d_{i,j}m_{j}f(r_{i,j}), \sigma^{2}),

where N() denotes the normal distribution, d_{i,j} is the tangential direction, i.e. the direction in which halo j bends the light of galaxy i, m_{j} is the mass of halo j, and f(r_{i,j}) is a decreasing function in the euclidean distance r_{i,j} between galaxy i and halo j.

After looking at the data I fixed the variance of the Gaussian distribution \sigma^{2} at 0.05. Like most competitors I also noticed that all skies seemed to have a single large halo, and that the other halos were much smaller. For the large halo I assigned the halo mass m a log-uniform distribution between 40 and 180, and I set f(r_{i,j}) = 1/\text{max}(r_{i,j},240). For the small halos I fixed the mass at 20, and I used f(r_{i,j}) = 1/\text{max}(r_{i,j},70). The resulting model is likely to be overly simplistic but it seems to capture most of the signal that is present in the data. In addition, keeping the model simple protected me against overfitting the data. Note that I assumed that the galaxy positions were independent of the halo positions, although it turns out this may not have been completely accurate.

After completing step 1 and 2, step 3 and 4 are simply a matter of implementation: I choose to use a simple random-walk Metropolis Hastings sampler to approximate the posterior distribution in step 3. The optimization in step 4 was done using standard gradient-based optimization, with random restarts to avoid local minima.

Like I remarked in the competition forums, the outcome of this competition was more noisy than is usual: final prediction accuracy was judged on a set of only 90 cases, with an evaluation metric that is very sensitive to small (angular) perturbations of the predictions. The public leaderboard standings were even more random, being based on only 30 cases. In fact, the 1.05 public score of my winning submission was only about average on the public leaderboard. All of this means I was very lucky indeed to win this competition. Nevertheless, runner-up Iain Murray seems to have taken a very similar approach, suggesting there is at least something to be said for looking at this kind of problem from a Bayesian perspective.

Finally, I would like to thank the organizers Dave and Tom, and sponsor Winton Capital for organizing a great competition. Looking at a problem different from the standard regression/classification problems was very refreshing.

I will post the Matlab code for my solution sometime later this week