Gaussian Symplectic Ensemble

In this post I define the Gaussian Symplectic Ensemble (GSE). Often the GSE is defined using quaternions; here I only use complex numbers. I give an accurate definition of the symplectic group, of self-dual matrices, of the GSE and give an algorithm to draw random matrices from the GSE.
The complex symplectic group is the set $$\begin{equation}\label{eq:20151211a} Sp(2 n , \mathbb{C} ) = \left\{ S \in M_{2 n} ( \mathbb{C} ) \Big| S^T J S = J\right\} \end{equation}$$ Hereby is $M_{2 n} ( \mathbb{C} )$ the set of complex matrices of size $2 n \times 2 n$ and $J$ is the $2 n \times 2 n$ matrix $$\begin{equation}\label{eq:20151211b} J = \left( \begin{array}{cc} 0& 1_n \\ -1_n & 0 \end{array} \right) \end{equation}$$ with $1_n$ the identity matrix of size $n \times n$. Notice that $Sp(2 n , \mathbb{C} )$ is defined with the transpose $S^T$, not with the Hermitian conjugate $S^\dagger$. One can verify that $Sp(2 n , \mathbb{C} )$ is a group with the usual matrix multiplication.
Definition of dual of a matrix
If $A$ is a complex matrix of size $2 n \times 2 n$, then its dual $A^D$ is defined as $$\begin{equation}\label{eq:20151211c} A^D = J^T A^T J \end{equation}$$ Notice again that this definition uses the transpose, not Hermitian conjugation. The following formulas are easily proved $$\begin{align*} (AB)^D &= B^D A^D\\ A^{DD} & = A \end{align*}$$ It is also easy to see that $$\begin{equation*} Sp(2 n , \mathbb{C} ) = \left\{ S \in M_{2 n} ( \mathbb{C} ) \Big| S^D S= 1 \right\} \end{equation*}$$ A self-dual matrix is a complex matrix $A$ of size $2 n \times 2 n$ which satisfies $$\begin{equation}\label{eq:20151211d} A^D = A \end{equation}$$
Definition of GSE
The GSE is the set of Hermitian, self-dual complex matrices $H$ of size $2 n \times 2 n$ with probability density function $$\begin{equation}\label{eq:20151211e} P(H) dH = \exp\left( -\frac{1}{2} \mathrm{tr} (H^2 ) \right) dH \end{equation}$$ where $dH$ is the product over all independent components $dH_{ij}$. The GSE is invariant under the action $$\begin{equation}\label{eq:20151211f} H \mapsto S^D H S \end{equation}$$ with $S \in USp(2 n ) = Sp(2 n , \mathbb{C} ) \cap U(2 n )$ . Hereby is $U(2 n )$ the group of unitary matrices of size $2 n \times 2 n$. It is easy to see that:

• If $H$ is self-dual, then $S^D H S$ is also self-dual.
• Because $S^DS = 1 = S^\dagger S$, it follows that $S^D = S^\dagger$. If $H$ is Hermitian, it thus follows that $S^D H S = S^\dagger H S$ is also Hermitian.

The action $\eqref{eq:20151211f}$ is therefore well defined. The group $USp(2 n )$ is called the unitary symplectic group.
Block form of $H$
If I write $H$ as a block matrix with $n \times n$ blocks, then $H$ is Hermitian and self-dual if and only if $H$ has the form $$\begin{equation}\label{eq:20151211g} H = \left( \begin{array}{cc} A& B \\ B^\dagger & A^T \end{array} \right) \end{equation}$$ with $A$ a Hermitian matrix of size $n \times n$ and $B$ a complex anti-symmetric matrix of size $n \times n$. The probability density function $\eqref{eq:20151211e}$ is Gaussian, with covariance $$\begin{equation}\label{eq:20151211h} \mathbb{E} \left[ H_{ij} H_{kl} \right] = \frac{1}{2} \left( \delta_{il} \delta_{jk} + J_{ik} J_{jl} \right) \end{equation}$$ It is easy to prove $\eqref{eq:20151211h}$ if one analyzes the blocks in $\eqref{eq:20151211g}$ case by case.
Algorithm to generate random matrices from the GSE
One can draw a random matrix from GSE as follows:

1. Draw a random matrix $A$ from the Gaussian Unitary Ensemble, with normalization $$\begin{equation*} \mathbb{E} \left[ A_{ij} A_{kl} \right] = \delta_{il} \delta_{jk} \end{equation*}$$
2. Calculate $H = \frac{1}{2} \left( A + A^D \right)$.

$H$ will then satisfy $\eqref{eq:20151211h}$