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

TL;DR

This paper presents a practical numerical scheme for solving Dyson Brownian motion to sample invariant matrix ensembles, demonstrating stability, fast convergence, and efficient computation methods.

Contribution

It introduces an explicit tamed Euler scheme for Dyson Brownian motion that effectively samples equilibrium measures despite singular Coulomb interactions.

Findings

Scheme is stable for time steps of O(1/N^2).

Convergence to equilibrium is exponential with rate O(1).

Sampling cost is optimized with fast multipole methods.

Abstract

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.