  We introduce a numerical algorithm to stochastically sample the dual fermion
perturbation series around the dynamical mean field theory, generating all
topologies of two-particle interaction vertices. We show results in the weak
and strong coupling regime of the half-filled Hubbard model in two dimensions,
illustrating that the method converges quickly where dynamical mean field
theory is a good approximation, and show that corrections are large in the
strong correlation regime at intermediate interaction. The fast convergence of
dual corrections to dynamical mean field results illustrates the power of the
approach and opens a practical avenue towards the systematic inclusion of
non-local correlations in correlated materials simulations. An analysis of the
frequency scale shows that only low-frequency propagators contribute
substantially to the diagrams, putting the inclusion of higher order vertices
within reach.
