Fast sampling from $\beta$-ensembles

Abstract

We investigate sampling $\beta$-ensembles with time complexity less than cubic in the cardinality of the ensemble. Following Dumitriu and Edelman (J Math Phys 43(11):5830–5847, 2002), we see the ensemble as the eigenvalues of a random tridiagonal matrix, namely a random Jacobi matrix. First, we provide a unifying and elementary treatment of the tridiagonal models associated with the three classical Hermite, Laguerre, and Jacobi ensembles. For this purpose, we use simple changes of variables between successive reparametrizations of the coefficients defining the tridiagonal matrix. Second, we derive an approximate sampler for the simulation of more general $\beta$-ensembles and illustrate how fast it can be for polynomial potentials. This method combines a Gibbs sampler on Jacobi matrices and the diagonalization of these matrices. In practice, even for large ensembles, only a few Gibbs passes suffice for the marginal distribution of the eigenvalues to fit the expected theoretical distribution. When the conditionals in the Gibbs sampler can be simulated exactly, the same fast empirical convergence is observed for the fluctuations of the largest eigenvalue. Our experimental results support a conjecture by Krishnapur et al. (Commun Pure Appl Math 69(1): 145–199, 2016), that the Gibbs chain on Jacobi matrices of size N mixes in $\mathcal {O}(\log N)$.

Appendices

Appendices

1.1 Appendix A

Proof

(of Lemma 1) For $n \in \mathbb {N}$, we note

$$\begin{aligned} \varDelta _n\left( x_{1}, \dots , x_{N} \right) \triangleq \begin{bmatrix} 1 &{} \cdots &{} 1 \\ x_1 &{} \cdots &{} x_N \\ &{} \vdots &{} \\ x_1^{n-1} &{} \cdots &{} x_N^{n-1} \end{bmatrix}. \end{aligned}$$

(44)

Then, for any $1\le n\le N$, the Cauchy–Binet formula (Lacanster and Tismenetsky 1985, Section 2.5) yields

$$\begin{aligned}&\det { \underline{H}_{2n-2} } = \det \left[ \sum _{k=1}^N w_k x_k^{i+j} \right] _{i,j=0}^{n-1} \nonumber \\&\quad = \det \left( \varDelta _n\left( x_{1}, \dots , x_{N} \right) \begin{bmatrix} w_1 &{} &{} \\ &{} \ddots &{} \\ &{} &{} w_N \end{bmatrix} \varDelta _n\left( x_{1}, \dots , x_{N} \right) ^{\textsf {T}} \right) \nonumber \\&\quad = \sum _{\left\{ i_1,\dots ,i_n \right\} \subset [N]} \left( \det \varDelta _n\left( x_{i_1}, \dots , x_{i_n} \right) \right) ^2 \prod _{k=1}^n w_{i_k} > 0. \end{aligned}$$

(45)

The particular case $n=N$ yields

$$\begin{aligned} \det { \underline{H}_{2N-2} } = \left( \det \varDelta _N\left( x_{1}, \dots , x_{N} \right) \right) ^2 \prod _{n=1}^N w_{n}. \end{aligned}$$

For $n>N$, (45) clearly shows that $\underline{H}_{2n-2}$ is rank deficient. The two other determinants are obtained similarly. $\square$

1.2 Appendix B

Proof

(of Proposition 4) Moments define monic OPs; see Proposition 3. By Favard’s theorem, here Theorem 5, monic OPs in turn define the atoms and weights of $\mu$ uniquely. Thus, $\phi$ is injective. Moreover, $\mathbb {R}_>^N\times S_N\subset \mathbb {R}^{2N-1}$ is open, and $\phi$ is $C^1$. By the classical inverse function theorem, see, e.g., Cartan (1971, Corollary 4.2.2), it is thus enough to show that the Jacobian of $\phi$ never vanishes.

To compute the partial derivative w.r.t. the nodes and weights, we write i-th moment of $\mu$ in two forms

$$\begin{aligned} m_i = \sum _{j= 1}^{N} w_j x_j^i = \sum _{j= 1}^{N- 1} w_j \left( x_j^i - x_N^i \right) + x_N^i. \end{aligned}$$

(46)

Thus, the Jacobian of $\phi$ reads

$$\begin{aligned}&\left| \frac{ \partial m_{1:2N- 1} }{ \partial x_{1:N}, w_{1:N- 1} } \right| = \left| \left[ \left[ \frac{\partial m_i}{\partial x_j} ~ \frac{\partial m_i}{\partial w_j} \right] _{j= 1}^{N- 1} ~ \left[ \frac{\partial m_i}{\partial x_N} \right] \right] _{i= 1}^{2N- 1} \right| \nonumber \\&\quad = \left| \left[ \left[ i w_j x_j^{i-1} ~~ x_j^{i} - x_N^{i} \right] _{j= 1}^{N- 1} ~ i w_N x_N^{i-1} \right] _{i= 1}^{2N- 1} \right| \nonumber \\&\quad = \left| \left[ \left[ (i-1)x_j^{i-2} ~~ x_j^{i-1} - x_N^{i-1} \right] _{j= 1}^{N- 1} ~ (i-1)x_N^{i-2} ~~ x_N^{i-1} \right] _{i= 1}^{2N} \right| \nonumber \\&\qquad \times \prod _{n=1}^N w_n \nonumber \\&\quad = \left| \left[ (i-1)x_j^{i-2} ~~ x_j^{i-1} \right] _{i=1,j= 1}^{2N,N} \right| \prod _{n=1}^N w_n. \end{aligned}$$

(47)

The last equality is obtained by adding the last column to all other even columns. The determinant in (47) is called a confluent Vandermonde determinant and has closed form expression given by $\prod _{i<j} (x_j-x_i)^{2\times 2} = \left( \det \varDelta _N(x_1,\dots ,x_N) \right) ^4$, see, e.g., Ha and Gibson (1980, Corollary 1, with $\eta _i \equiv 2$). In particular, (47) never vanishes on $\mathbb {R}_>^N\times S_N$. $\square$

1.3 Appendix C

Proof

(of Proposition 5) Using Theorem 5 and Proposition 4, $\rho = \psi \circ \phi ^{-1}$, so that $\rho$ is bijective. As in the proof of Proposition 4, we apply the inverse function theorem (Cartan 1971, Corollary 4.2.2), but this time to $\rho ^{-1}$. We first note that $\mathbb {R}^N \times (0, +\infty )^{N-1}\subset \mathbb {R}^{2N-1}$ is open. It is thus enough to show that $\rho ^{-1}$ is $C^1$ and that its Jacobian never vanishes. To this end, we borrow an elegant lattice path representation of the recurrence relations for OPs from Hardy (2017, Equation 1.8). This allows us to express the successive moments as polynomials in the recurrence coefficients.

To provide intuition, we first compute the first few moments by hand, recursively applying the recurrence relation (12). It comes

(48)

Fig. 7

The lattice path of Hardy (2017) used to compute $m_n=\left\langle x^nP_0, P_0 \right\rangle$ is displayed in (a). The paths used for the computation of $m_3$ (48) are highlighted in b–e with the corresponding weight as caption

More generally, when computing $m_k = \left\langle x^k P_0, P_0 \right\rangle$, successive applications of the recurrence relation (12) allow to decrease the power of x from k to 0 until each term in the expansion is proportional to the inner product of $P_0=1$ with another monic OP. The only nonzero such inner product is $\left\langle P_0,P_0 \right\rangle =1$. Consequently, each nonzero term in the final expansion of $m_k$ corresponds to a path of length at most k that leaves from the lower left corner of the graph in Fig. 7a and ends up on the bottom row. In between, the path has to remain above the bottom row, and can only move North-East, East, or South-East. Each edge corresponds to picking one of the three terms in the recurrence relation (12). For example, the expansion of $m_3$ in (48) corresponds to three such paths, shown in orange in Fig. 7. The product of the coefficients along each path forms the resulting term in the expansion.

In the end, odd moments $m_{2i-1}$, resp. even moments $m_{2i}$, are the sum of the weights of the paths below the i-th red, respectively blue path, counting from the bottom left. More precisely,

$$\begin{aligned} m_{2i-1}&= { a_{i} \prod _{k= 1}^{i-1} b_{k} } + f_1( a_{1:i-1}, b_{1:i-2}), \text { and} \nonumber \\ m_{2i}&= { \prod _{k= 1}^{i} b_{k} } + f_2( a_{1:i}, b_{1:i-1}). \end{aligned}$$

(49)

Thus, the Jacobian matrix is triangular with determinant

$$\begin{aligned} \left| \frac{ \partial m_{1:2N-1} }{ \partial a_{1:N}, b_{1:N-1} } \right| =\prod _{i= 1}^N \frac{\partial m_{2i-1}}{\partial a_{i}} \prod _{i= 1}^{N-1} \frac{\partial m_{2i}}{\partial b_{i}} ,\end{aligned}$$

and the formulation (49) yields

$$\begin{aligned}&\left| \frac{ \partial m_{1:2N-1} }{ \partial a_{1:N}, b_{1:N-1} } \right| = \prod _{i= 1}^N \prod _{k= 1}^{i-1} b_{k} \prod _{i= 1}^{N-1} \prod _{k= 1}^{i-1} b_{k}\\&\quad = \frac{ \left[ \prod _{i= 1}^{N} \prod _{k= 1}^{i-1} b_{k} \right] ^2 }{ \prod _{k= 1}^{N-1} b_{k} } = \frac{ \left[ \prod _{n= 1}^{N-1} b_{n}^{N- n} \right] ^2 }{ \prod _{n= 1}^{N-1} b_{n} },\end{aligned}$$

which does not vanish since all $b_n$s are positive by construction. Finally, the last equality in (23) follows from Lemma 2. $\square$

1.4 Appendix D

Proof

(of Lemma 6) First, combine the result of Lemma 4 and the definition of the canonical moments in (35) to get

$$\begin{aligned} \prod _{n=1}^N x_n&= \xi _1 \prod _{n=2}^N \xi _{2n-1} \\&= c_1 \prod _{n=2}^N (1-c_{2n-2})c_{2n-1} = \prod _{n=1}^N c_{2n-1} (1-c_{2n}). \end{aligned}$$

Then, Lemma 1 yields

$$\begin{aligned} \prod _{n= 1}^N (1-x_n) = \frac{\det { \overline{H}_{2N-1} } }{\det { \underline{H}_{2N-2} } } \cdot \end{aligned}$$

(50)

The denominator can be expressed using the $\xi _{1:2N-1}$ parameters,

(51)

For the numerator, we follow Dette and Studden (1997, Theorem 1.4.10) who introduced additional quantities $\gamma _{1:2N-1}$ to get

$$\begin{aligned} \det { \overline{H}_{2N-1} } = \gamma _{1}^N \prod _{n=1}^{N- 1} \left[ \gamma _{2n}\gamma _{2n+1} \right] ^{N-n}, \end{aligned}$$

(52)

where $\xi _1 = c_1, \gamma _1 = 1 - c_1$, and for any $1\le n \le 2N-1$,

$$\begin{aligned} \xi _n = (1-c_{n-1})c_n, \text { and } \gamma _{n} = c_{n-1} (1-c_n). \end{aligned}$$

We plug these results back into (50) and conclude that

$\square$