Random Thoughts in a Random Universe

Convergence of the Affine Invariant Ensemble Sampler

The Affine-Invariant Ensemble Sampler (AIES) is a Monte-Carlo method particularly popular in astrophysics. The paper describing its widely used Python implementation emcee is highly cited and the method is particularly popular because it has — aside from the number of walkers, i.e., the number of chains — no tunable parameters.

Huijser, Goodman, and Brewer recently highlighted some convergence problems in AIES. This is illustrated in the case of correlated Gaussians, specifically in the auto-regressive process of order 1, aka AR(1). This blog post sets out to test whether a simple convergence criterion can discover these problems.

The Autoregressive (1) Process

Huijser, Goodman, and Brewer use the so-called autoregressive process of order 1 (aka AR(1)) to highlight slow convergence of the AIES algorithm when the dimensionality of the problem exceeds a few and correlations between parameters exist. We start by reproducing their findings.

In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pylab as plt

import emcee
import corner

We set up the covariance matrix for the AR(1) process in 10, 50, and 100 dimensions.

In [2]:
ndim_arr = [10, 50, 100]
α = 0.9
cov_arr = []
for ndim in ndim_arr:
    idx = np.repeat(np.arange(ndim), ndim).reshape(ndim, ndim)
    Σ = α**abs(idx - idx.T)  
    cov_arr.append(Σ)    

Let’s first illustrate the correlated distribution for the 10-dimensional case.

In [3]:
X = np.random.multivariate_normal(np.zeros(ndim_arr[0]), cov_arr[0], 50000)
figure = corner.corner(X)