Numerical solution of Dyson Brownian motion and a sampling scheme for invariant matrix ensembles

The Dyson Brownian Motion (DBM) describes the stochastic evolution of $N$ points on the line driven by an applied potential, a Coulombic repulsion and identical, independent Brownian forcing at each point. We use an explicit tamed Euler scheme to numerically solve the Dyson Brownian motion and sample the equilibrium measure for non-quadratic potentials. The Coulomb repulsion is too singular for the SDE to satisfy the hypotheses of rigorous convergence proofs for tamed Euler schemes. Nevertheless, in practice the scheme is observed to be stable for time steps of O(1/N^2) and to relax exponentially fast to the equilibrium measure with a rate constant of O(1) independent of $N$. Further, this convergence rate appears to improve with $N$ in accordance with O(1/N) relaxation of local statistics of the Dyson Brownian motion. This allows us to use the Dyson Brownian motion to sample $N\times N$ Hermitian matrices from the invariant ensembles. The computational cost of generating $M$ independent samples is O(MN^4) with a naive scheme, and O(MN^3\log N) when a fast multipole method is used to evaluate the Coulomb interaction.