<jupyter_start><jupyter_text>$\beta$-EnsemblesA $\beta$-Ensemble is a (unordered) set of points with joint probability distribution $$(x_1,\dots,x_N) \sim     \frac{1}{Z_{N,\beta}}	\left|\Delta(x_1,\dots,x_N)\right|^{\beta}	\prod_{i= 1}^N 		 \mu_{\text{ref}}(\text{d}x_i)$$We consider the following reference measures $\mu_{\text{ref}}$- Normal distribution $\mathcal{N}(\mu, \sigma)$, associated to the Hermite ensemble.$$    \frac{1}{\sqrt{2\pi\sigma^2}}        e^{- \frac{1}{2\sigma^2} (x-\mu)^2}    1_{\mathbb{R}}(x)    dx$$- Gamma distribution $\Gamma(k, \theta)$, associated to the Laguerre ensemble.$$ \frac{1}{\theta^k\Gamma(k)}      x^{k-1}      e^{- \frac{x}{\theta}}     1_{\mathbb{R}^+}(x)    dx$$- Beta distribution $\operatorname{Beta}(a,b)$, associated to the Jacobi ensemble.$$\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}        x^{a-1}        (1-x)^{b-1}    1_{(0,1)}(x)    dx$$- Uniform distribution on the unit circle $\mathcal{U}_{\mathbb{S}^1}$, associated to the Circular ensemble.$$\frac{1}{2\pi}      1_{|x|=1}    d x$$- Complex Gaussian distribution, associated to the Ginbre ensemble ($\beta=2$).$$e^{ - \frac{1}{2}|z|^2 }d z$$<jupyter_code>%pylab inline

%load_ext autoreload
%autoreload 2

import os
import sys
sys.path.insert(0, os.path.abspath('../dppy'))

from beta_ensembles import *
supported_ensembles = ('hermite', 'laguerre', 'jacobi', 'circular', 'ginibre')
err_print = ['Supported OPEs:']
err_print.extend(['- {}'.format(OPE) for OPE in supported_ensembles])

'\n'.join(err_print)<jupyter_output><empty_output><jupyter_text>Ensembles on the real lineRandom matrix theory...Some matrix models exist for $\beta=1,2,4$ casesTridiagonal models for Hermite and Laguerre[<NAME> (2002)](https://sites.math.washington.edu/~dumitriu/JMathPhys_43_5830.pdf)  Hermite$$(x_1,\dots,x_N) \sim 	\left|\Delta(x_1,\dots,x_N)\right|^{\beta}	\prod_{i= 1}^N 		e^{- \frac{1}{2}\frac{x_i^2}{2}} 		% \indic_{\bbR}(x_i)	\ d x_i$$ Full matrix model GOE ($\beta=1$)$$	X \sim \mathcal{N}_{N, N}(0,1)	\qquad	A = \frac{X+X^{\top}}{\sqrt{2}}$$ GUE  ($\beta=2$)$$	X \sim \mathcal{N}_{N, N}(0,1) + i~ \mathcal{N}_{N, N}(0,1)	\qquad	A = \frac{X+X^{\dagger}}{\sqrt{2}}$$ GSE  ($\beta=4$)$$	\begin{cases}        X \sim \mathcal{N}_{N, M}(0,1) + i~ \mathcal{N}_{N, M}(0,1)\\        Y \sim \mathcal{N}_{N, M}(0,1) + i~ \mathcal{N}_{N, M}(0,1)    \end{cases}    \qquad	A = \begin{bmatrix}            X    & Y   \\            -Y^* & X^*        \end{bmatrix}    \quad A = \frac{X+X^{\dagger}}{\sqrt{2}}$$Normalization $\sqrt{\beta N}$ to concentrate as the semi-circle law.<jupyter_code>ensemble_name, beta = 'hermite', 2
hermite = BetaEnsemble(ensemble_name, beta=beta)
print(hermite) # or hermite.info()
hermite.flush_samples()
mode = 'full'

# Plot
hermite_params = {'N':30}
hermite.sample(mode, **hermite_params)
hermite.plot(normalization=True)

# Histo
hermite_params = {'N':2000}
hermite.sample(mode, **hermite_params)
hermite.hist(normalization=True)<jupyter_output><empty_output><jupyter_text>Tridiagonal model ($\beta>0$)    [<NAME> (2002), II-C](https://sites.math.washington.edu/~dumitriu/JMathPhys_43_5830.pdf)$$\begin{bmatrix}    \alpha_1 		& \sqrt{\beta_2}& 0         	&    0      & 0       	\\    \sqrt{\beta_2}  & \alpha_2 		& \sqrt{\beta_3}&    0      & 0       	\\        0    		& \ddots   		& \ddots      	& \ddots	& 0       	\\        0 			& 0 			& \sqrt{\beta_{N-1}} & \alpha_{N- 1}   & \sqrt{\beta_{N}} \\        0   		& 0        & 0           	& \sqrt{\beta_N} 	& \alpha_{N}\end{bmatrix}$$$$\alpha_{i}    \sim \mathcal{N}(0, \sqrt{2})    \quad\text{and}    \quad\beta_{i+1}    \sim \chi_{\beta(N - i)}^2$$ More generally for a $\mathcal{N}(\mu, \sigma)$ reference measure$$(x_1,\dots,x_N) \sim 	\left|\Delta(x_1,\dots,x_N)\right|^{\beta}	\prod_{i= 1}^N 		e^{- \frac{(x_i-\mu)^2}{2\sigma^2}} 		% \indic_{\bbR}(x_i)	\ d x_i$$$$\alpha_{i}    \sim \mathcal{N}(\mu,\sigma)    \quad\text{and}    \quad\beta_{i+1}    \sim \Gamma\left(\frac{\beta}{2}(N - i), \sigma^2\right)$$Recall $\Gamma(\frac{k}{2}, 2)\equiv \chi_k^2$<jupyter_code>ensemble_name, beta = 'hermite', 4.1
hermite = BetaEnsemble(ensemble_name, beta=beta)
print(hermite) # or hermite.info()
hermite.flush_samples()
mode = 'banded'

# Plot
hermite_params = {'loc':0.0, 'scale':1.0, 'size':30}
# To match the full matrix model use
# loc = 0.0
# scale = np.sqrt(2)
hermite.sample(mode, **hermite_params)
hermite.plot(normalization=True)

# Histo
hermite_params = {'loc':0.0, 'scale':1.0, 'size':2000}
hermite.sample(mode, **hermite_params)
hermite.hist(normalization=True)<jupyter_output><empty_output><jupyter_text>Laguerre[<NAME> (2002)](https://sites.math.washington.edu/~dumitriu/JMathPhys_43_5830.pdf)$$(x_1,\dots,x_N) \sim 	\left|\Delta(x_1,\dots,x_N)\right|^{\beta}    %	\prod_{i= 1}^N         x_i^{\frac{\beta}{2}(M-N+1)-1}		e^{- \frac12 x_i}		% \indic_{\bbR}(x_i)	\ d x_i$$ Full matrix model LOE ($\beta=1$)$$	X \sim \mathcal{N}_{N, M}(0,1)	\qquad	A = XX^{\top}$$ LUE  ($\beta=2$)$$	X \sim \mathcal{N}_{N, M}(0,1) + i~ \mathcal{N}_{N, M}(0,1)	\qquad	A = XX^{\dagger}$$ LSE  ($\beta=4$)$$	\begin{cases}        X \sim \mathcal{N}_{N, M}(0,1) + i~ \mathcal{N}_{N, M}(0,1)\\        Y \sim \mathcal{N}_{N, M}(0,1) + i~ \mathcal{N}_{N, M}(0,1)    \end{cases}    \qquad	A = \begin{bmatrix}            X    & Y   \\            -Y^* & X^*        \end{bmatrix}    \quad A = A A^{\dagger}$$Normalization $\beta M$ to concentrate as Marcenko-Pastur law<jupyter_code>ensemble_name, beta = 'laguerre', 4
laguerre = BetaEnsemble(ensemble_name, beta=beta)
print(laguerre) # or laguerre.info()
laguerre.flush_samples()
mode='full'

# Plot
laguerre_params = {'M':40, 'N':30}
laguerre.sample(mode, **laguerre_params)
laguerre.plot(normalization=True)

# Histo
laguerre_params = {'M':2500, 'N':2000}
laguerre.sample(mode, **laguerre_params)
laguerre.hist(normalization=True)<jupyter_output><empty_output><jupyter_text>Tridiagonal model ($\beta>0$)[<NAME> (2002), III-B](https://sites.math.washington.edu/~dumitriu/JMathPhys_43_5830.pdf)$$\begin{bmatrix}    \alpha_1 		& \sqrt{\beta_2}& 0         	&    0      & 0       	\\    \sqrt{\beta_2}  & \alpha_2 		& \sqrt{\beta_3}&    0      & 0       	\\        0    		& \ddots   		& \ddots      	& \ddots	& 0       	\\        0 			& 0 			& \sqrt{\beta_{N-1}} & \alpha_{N- 1}   & \sqrt{\beta_{N}} \\        0   		& 0        & 0           	& \sqrt{\beta_N} 	& \alpha_{N}\end{bmatrix}=\begin{bmatrix}    \sqrt{\xi_1} 	& 	 			&          		&        	\\    \sqrt{\xi_2}  	& \sqrt{\xi_3}	&      			&        	\\         			& \ddots		& \ddots 		&   		\\          	 		&       		& \sqrt{\xi_{2N-2}} & \sqrt{\xi_{2N-1}}\end{bmatrix}\begin{bmatrix}    \sqrt{\xi_1} 	& \sqrt{\xi_2}	&          	&        					 \\      				& \sqrt{\xi_3}	& \ddots  	&        					 \\         			& 				& \ddots 	& \sqrt{\xi_{2N-2}} \\          	 		&       		&			& \sqrt{\xi_{2N-1}}\end{bmatrix}$$$$\begin{aligned}    \alpha_1 &= x_1 	                \quad & \\    \alpha_k &= \xi_{2k-2}+\xi_{2k-1}    \quad &\beta_{k+1} &= \xi_{2k-1}\xi_{2k}\end{aligned}$$$$\xi_{2i-1}    \sim \chi_{\beta(M - i + 1)}^2    \quad\text{and}    \quad\xi_{2i}    \sim \chi_{\beta(N - i)}^2$$ More generaly for a $\Gamma(k,\theta)$ reference measure$$(x_1,\dots,x_N) \sim 	\left|\Delta(x_1,\dots,x_N)\right|^{\beta}    %	\prod_{i= 1}^N         x_i^{k-1}		e^{- \frac{x_i}{\theta}} 		% \indic_{\bbR}(x_i)	\ d x_i$$$$\xi_{2i-1}    \sim \Gamma\left(\frac{\beta}{2}(N - i) + k, \theta \right)    \quad\text{and}    \quad\xi_{2i}    \sim \Gamma\left(\frac{\beta}{2}(N - i), \theta \right)$$Recall $\Gamma(\frac{k}{2}, 2)\equiv \chi_k^2$<jupyter_code>ensemble_name, beta = 'laguerre', 4
laguerre = BetaEnsemble(ensemble_name, beta=beta)
print(laguerre) # or laguerre.info()
laguerre.flush_samples()
mode='banded'

# Plot
laguerre_params = {'shape':22, 'scale':2.0, 'size':30}
# To match the full matrix model use
# shape = 0.5*beta*(M-N+1)
# scale = 2.0
laguerre.sample(mode, **laguerre_params)
laguerre.plot(normalization=True)

# Histo
laguerre_params = {'shape':10000, 'scale':2.0, 'size':2000}
laguerre.sample(mode, **laguerre_params)
laguerre.hist(normalization=True)<jupyter_output><empty_output><jupyter_text>Jacobi $$(x_1,\dots,x_N) \sim 	\left|\Delta(x_1,\dots,x_N)\right|^{\beta}    %	\prod_{i= 1}^N         x_i^{\frac{\beta}{2}(M_1-N+1)-1}        (1-x_i)^{\frac{\beta}{2}(M_2-N+1)-1}		% \indic_{\bbR}(x_i)	\ d x_i$$Concentrates as Wachter law Full matrix model JOE ($\beta=1$)$$	\begin{cases}		X \sim \mathcal{N}_{N, M_1}(0,1)\\		Y \sim \mathcal{N}_{N, M_2}(0,1)\\	\end{cases}	\qquad	A = XX^{\top}\left(XX^{\top} + YY^{\top}\right)^{-1}$$ JUE ($\beta=2$)$$	\begin{cases}		X \sim \mathcal{N}_{N, M_1}(0,1) + i~ \mathcal{N}_{N, M_1}(0,1)\\		Y \sim \mathcal{N}_{N, M_2}(0,1) + i~ \mathcal{N}_{N, M_2}(0,1)\\	\end{cases}	\qquad	A = XX^{\dagger}\left(XX^{\dagger} + YY^{\dagger}\right)^{-1}$$ JSE ($\beta=4$)$$	\begin{cases}		X_1 \sim \mathcal{N}_{N, M_1}(0,1) + i~ \mathcal{N}_{N, M_1}(0,1)\\        X_2 \sim \mathcal{N}_{N, M_1}(0,1) + i~ \mathcal{N}_{N, M_1}(0,1)\\        Y_1 \sim \mathcal{N}_{N, M_2}(0,1) + i~ \mathcal{N}_{N, M_2}(0,1)\\		Y_2 \sim \mathcal{N}_{N, M_2}(0,1) + i~ \mathcal{N}_{N, M_2}(0,1)	\end{cases}	\qquad    \begin{cases}        X = \begin{bmatrix}                X_1    & X_2   \\                -X_2^* & X_1^*            \end{bmatrix}\\        Y = \begin{bmatrix}                Y_1    & Y_2   \\                -Y_2^* & Y_1^*            \end{bmatrix}    \end{cases}    \qquad	A = XX^{\dagger}\left(XX^{\dagger} + YY^{\dagger}\right)^{-1}$$<jupyter_code>ensemble_name, beta = 'jacobi', 2
jacobi = BetaEnsemble(ensemble_name, beta=beta)
print(jacobi) # or jacobi.info()
jacobi.flush_samples()
mode='full'

# Plot
jacobi_params = {'M_1':60, 'M_2':40, 'N':20}
jacobi.sample(mode, **jacobi_params)
jacobi.plot(normalization=True)

# Histo, see convergence towards the arcsine distribution
jacobi_params = {'M_1':2500, 'M_2':2200, 'N':2000}
jacobi.sample(mode, **jacobi_params)
jacobi.hist(normalization=True)<jupyter_output><empty_output><jupyter_text>Tridiagonal model ($\beta>0$)see [(Killip-Nenciu, Theroem 2)](https://arxiv.org/pdf/math/0410034.pdf)$$\begin{bmatrix}    \alpha_1 		& \sqrt{\beta_2}& 0         	&    0      & 0       	\\    \sqrt{\beta_2}  & \alpha_2 		& \sqrt{\beta_3}&    0      & 0       	\\        0    		& \ddots   		& \ddots      	& \ddots	& 0       	\\        0 			& 0 			& \sqrt{\beta_{N-1}} & \alpha_{N- 1}   & \sqrt{\beta_{N}} \\        0   		& 0        & 0           	& \sqrt{\beta_N} 	& \alpha_{N}\end{bmatrix}$$$$\begin{aligned}    \alpha_1 &= \xi_1 	                \quad & \\    \alpha_k &= \xi_{2k-2}+\xi_{2k-1}    \quad &\beta_{k+1} &= \xi_{2k-1}\xi_{2k}\end{aligned}$$$$\begin{aligned}    \xi_1 &= c_1 	           \quad &\gamma_1 &= 1-c_1 \\    \xi_k &= (1-c_{k-1})c_k    \quad &\gamma_k &= c_{k-1}(1-c_k)\end{aligned}$$$$(x_1,\dots,x_N) \sim 	\left|\Delta(x_1,\dots,x_N)\right|^{\beta}    %	\prod_{i= 1}^N         x_i^{\frac{\beta}{2}(M_1-N+1)-1}        (1-x_i)^{\frac{\beta}{2}(M_2-N+1)-1}		% \indic_{\bbR}(x_i)	\ d x_i$$$$c_{2i-1}    \sim     \operatorname{Beta}    \left(         \frac{\beta}{2}(M_1-i+1),         \frac{\beta}{2}(M_2-i+1)     \right)    \quad\text{and}    \quadc_{2i}    \sim    \operatorname{Beta}    \left(         \frac{\beta}{2} (N-i),         \frac{\beta}{2} (M_1+M_2-N-i+1)     \right)$$ More generally for a $\operatorname{Beta}(a,b)$ reference measure$$(x_1,\dots,x_N) \sim 	\left|\Delta(x_1,\dots,x_N)\right|^{\beta}    %	\prod_{i= 1}^N         x_i^{a-1}        (1-x_i)^{b-1}		% \indic_{\bbR}(x_i)	\ d x_i$$$$c_{2i-1}    \sim     \operatorname{Beta}    \left(         \frac{\beta}{2}(N-i) + a,         \frac{\beta}{2}(N-i) + b     \right)    \quad\text{and}    \quadc_{2i}    \sim    \operatorname{Beta}    \left(         \frac{\beta}{2} (N-i),         \frac{\beta}{2} (N-i-1) + a + b     \right)$$<jupyter_code>ensemble_name, beta = 'jacobi', 4.5
jacobi = BetaEnsemble(ensemble_name, beta=beta)
print(jacobi) # or jacobi.info()
jacobi.flush_samples()
mode = 'banded'

# Plot
jacobi_params = {'a':30, 'b':10, 'size':20}
# To match the full matrix model use
# a, b = 0.5*beta*(M_1-N+1), 0.5*beta*(M_2-N+1)
jacobi.sample(mode, **jacobi_params)
jacobi.plot(normalization=True)

# Histo, see convergence towards the arcsine distribution
jacobi_params = {'a':300, 'b':100, 'size':300}
jacobi.sample(mode, **jacobi_params)
jacobi.hist(normalization=True)<jupyter_output><empty_output><jupyter_text>Classical Ensembles on the circle and in the plane  Circular ensemble$$\left|\Delta(e^{i \theta_1 },\dots, e^{i \theta_N})\right|^{\beta}    \prod_{j = 1}^N \frac{1}{2\pi} \mathbf{1}_{[0,2\pi]} (\theta_j) d\theta_j$$ Full matrix model COE ($\beta=1$)  Uniform measure i.e. Haar measure on orthogonal matrices $\mathbb{O}_N$:  $U^{\top}U = I_N$ Via QR algorithm, see [Mezzadri, Sec 5](https://arxiv.org/pdf/math-ph/0609050.pdf)```pythonA = np.random.randn(N, N)Q, R = la.qr(A)d = np.diagonal(R)U = np.multiply(Q, d/np.abs(d), Q)return eigvals(U)``` The Hermite way$	X \sim \mathcal{N}_{N, N}(0,1)\\	A = X+X^{\top}       = U^{\top}\Lambda U\\    eigvals(U)$ CUE  ($\beta=2$)  Uniform measure i.e. Haar measure on unitary matrices $\mathbb{U}_N$:   $U^{\dagger}U = I_N$ Via QR algorithm, see [Mezzadri, Sec 5](https://arxiv.org/pdf/math-ph/0609050.pdf)```pythonA = np.random.randn(N, N) + 1j*np.random.randn(N, N)A /= np.sqrt(2.0)Q, R = la.qr(A)d = np.diagonal(R)U = np.multiply(Q, d/np.abs(d), Q)return eigvals(U)``` The Hermite way$	X \sim \mathcal{N}_{N, N}(0,1) + i~ \mathcal{N}_{N, N}(0,1)\\    A = X+X^{\dagger}      = U^{\dagger}\Lambda U\\    eigvals(U)$ CSE  ($\beta=4$) Uniform measure i.e. Haar measure  onsymplectic matrices $\mathbb{U}\operatorname{Sp}_{2N}$:   $U^{\dagger}U = I_N$$	\begin{cases}        X \sim \mathcal{N}_{N, M}(0,1) + i~ \mathcal{N}_{N, M}(0,1)\\        Y \sim \mathcal{N}_{N, M}(0,1) + i~ \mathcal{N}_{N, M}(0,1)    \end{cases}\\	A = \begin{bmatrix}            X    & Y   \\            -Y^* & X^*        \end{bmatrix}    \quad A = X+X^{\dagger}            = U^{\dagger} \Lambda U\\    eigvals(U)$<jupyter_code>ensemble_name, beta = 'circular', 4
circular = BetaEnsemble(ensemble_name, beta=beta)
circular.info()
circular.flush_samples()
mode='full'

# Plot
circular_params = {'N':20, 'haar_mode':'hermite'} # 'haar_mode':'QR'
circular.sample(mode, **circular_params)
circular.plot(normalization=True)

# Histo
circular_params = {'N':1000, 'haar_mode':'hermite'} # 'haar_mode':'QR'
circular.sample(mode, **circular_params)
circular.hist(normalization=True)<jupyter_output><empty_output><jupyter_text>Quindiagonal modelsee [(Killip-Nenciu, Theroem 1)](https://arxiv.org/pdf/math/0410034.pdf)$$\left|\Delta(e^{i \theta_1 },\dots, e^{i \theta_N})\right|^{\beta}    \prod_{j = 1}^N \frac{1}{2\pi} \mathbf{1}_{[0,2\pi]} (\theta_j) d\theta_j$$For $\nu\geq2$ an integer, if $v$ is chosen uniformly at random from the unit sphere $\mathbb{S}^{\nu} \in \mathbb{R}^{\nu+1}$, then $v_1 + iv_2\sim \Theta_{\nu}$Given $\beta\in \mathbb{N}^*$, let $\alpha_k\sim \Theta_{\beta(N-k-1)+1}$ be independent random variablesfor $0\leq k\leq N-1$, $\rho_k = \sqrt{1-|\alpha_k|^2}$, and define $$\Xi_k = \begin{bmatrix}\overline{\alpha}_k & \rho_k\\\rho_k   & -\alpha_k\end{bmatrix}, \quad 0\leq k\leq N-2, \quad \text{with} \quad \Xi_{-1} = [1]\quad \text{and} \quad\Xi_{N-1} = [\overline{\alpha}_{N-1}]$$$$L = \operatorname{diag}[\Xi_0,\Xi_2,\dots]\quad \text{and} \quadM = \operatorname{diag}[\Xi_{-1},\Xi_1,\Xi_3\dots]$$$LM$ and $ML$ have same eigenvalue distribution<jupyter_code>ensemble_name, beta = 'circular', 10
circular = BetaEnsemble(ensemble_name, beta=beta)
circular.info()
mode = 'banded'
circular.flush_samples()

# Plot
circular_params = {'size':20}
circular.sample(mode, **circular_params)
circular.plot(normalization=True)

# Histo
circular_params = {'size':1000}
circular.sample(mode, **circular_params)
circular.hist(normalization=True)<jupyter_output><empty_output><jupyter_text>Ginibre$$\left|\Delta(z_1,\dots,z_N)\right|^{2}\prod_{i = 1}^N e^{ - \frac{1}{2}|z_i|^2 }d z_i$$   Full matrix model$$A \sim \frac{1}{\sqrt{2}} \left( \mathcal{N}_{N,N}(0,1) + i~ \mathcal{N}_{N, N}(0,1) \right)$$Nomalization $\sqrt{N}$ to concentrate in the unit circle<jupyter_code>ensemble_name = 'ginibre' # No need for beta (=2)
ginibre = BetaEnsemble(ensemble_name)

ginibre.info()
ginibre.flush_samples()
ginibre_params = {'N':100, 'mode':'full'}
ginibre.sample(**ginibre_params)

ginibre.info()
ginibre.plot(normalization=True)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.22696969696969696, 'fraction_numerical': 0.03666666666666667, 'mean_word_length': 3.3027379400260757, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 48, 'https://': 8, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '33238', 'license': 'MIT', 'n_tokens_mistral': 7283, 'n_tokens_neox': 6891, 'n_words': 1184, 'path': 'notebooks/beta_ensembles.ipynb', 'repo_name': 'claunay/DPPy', 'size': '16544'}