# The conditional DPP approach to random matrix distributions

Alan Edelman  Department of Mathematics and Computer Science & AI Laboratory,
Massachusetts Institute of Technology, Cambridge, Massachusetts, 02139
<EMAIL_ADDRESS>and  Sungwoo Jeong  Department of Mathematics, Massachusetts
Institute of Technology, Cambridge, Massachusetts, 02139<EMAIL_ADDRESS>
###### Abstract.

We present the conditional determinantal point process (DPP) approach to
obtain new (mostly Fredholm determinantal) expressions for various eigenvalue
statistics in random matrix theory. It is well-known that many (especially
$\beta=2$) eigenvalue $n$-point correlation functions are given in terms of
$n\times n$ determinants, i.e., they are continuous DPPs. We exploit a derived
kernel of the conditional DPP which gives the $n$-point correlation function
conditioned on the event of some eigenvalues already existing at fixed
locations.

Using such kernels we obtain new determinantal expressions for the joint
densities of the $k$ largest eigenvalues, probability density functions of the
$k^{\text{th}}$ largest eigenvalue, density of the first eigenvalue spacing,
and more. Our formulae are highly amenable to numerical computations and we
provide various numerical experiments. Several numerical values that required
hours of computing time could now be computed in seconds with our expressions,
which proves the effectiveness of our approach.

We also demonstrate that our technique can be applied to an efficient sampling
of DR paths of the Aztec diamond domino tiling. Further extending the
conditional DPP sampling technique, we sample Airy processes from the extended
Airy kernel. Additionally we propose a sampling method for non-Hermitian
projection DPPs.

##  1\. Introduction

###  1.1. The Conditional DPP Method

This paper shows that conditional determinantal point processes (DPP) can be
exploited in novel ways to create interesting expressions and yield highly
efficient algorithms for computations in random matrix theory and beyond. We
call this the conditional DPP approach. Figure 1 presents a gallery of
examples created by these new expressions.

|
---|---
Section 3.4, Equation (1.5) | Section 3.4
|
Section 3.4.2 | Section 3.4.2, Equation (3.13)
|  Correlation coefficient (soft-edge)
$\rho(\lambda_{1},\lambda_{2})=0.50564723159...$
$\begin{array}[]{|c|}\hline\cr\text{Computation time}\\\ \hline\cr\text{16
hours (2010, \cite[cite]{[\@@bibref{}{bornemann2010numerical}{}{}]})}\\\
\downarrow\\\ \text{2 minutes (Section~{}\ref{sec:klargest}})\\\
\hline\cr\end{array}$
Sections 3.2, 3.3, Eq (1.6), (3.10) | Section 3.4.1
Figure 1. Plots and numerical results enabled by the conditional DPP approach.
Note that $\alpha$ is the parameter of the LUE.

###  1.2. Technical Background

Determinantal representations frequently arise in random matrix theory,
especially in the study of eigenvalues of $\beta=2$ (complex) random matrices.
A number of random matrix $n$-point eigenvalue correlation functions [8] are
given in terms of the following determinantal formula

(1.1)
$p(x_{1},\dots,x_{n})=\det\left(\left[K(x_{i},x_{j})\right]_{i,j=1\dots,n}\right),$

for some corresponding kernel $K$. A basic example is the $N\times N$ Gaussian
unitary ensemble (GUE) with the Hermite kernel $K=K_{\text{Herm}}^{(N)}$
defined as

(1.2)
$K_{\text{Herm}}^{(N)}(x,y)=\sum_{i=0}^{N-1}\phi_{i}(x)\phi_{i}(y)=\sqrt{\frac{N}{2}}\frac{\phi_{N}(x)\phi_{N-1}(y)-\phi_{N-1}(x)\phi_{N}(y)}{x-y},$

with the Hermite functions
$\phi_{j}(x)=\exp(-x^{2}/2)H_{j}(x)/(2^{j}\sqrt{\pi}j!)^{1/2}$, where
$H_{j}$’s are the Hermite polynomials.

Defined by such determinantal representations, a point process is called a
determinantal point process, if a random element $\mathcal{J}$ satisfies for
any subset $S$ of the ground index (point) set $\mathcal{G}$,

(1.3)
$\mathbb{P}(\text{$S\subset\mathcal{J}$})=\det\left(\left[K(x,y)\right]_{x,y\in
S}\right),$

given a kernel (or a matrix, in discrete cases) $K$. As we mentioned above,
eigenvalue statistics of $\beta=2$ orthogonal polynomial ensembles (and their
limits) are typical examples of the continuous DPP [13, 25]. In this paper the
kernel $K$ is not restricted to symmetric or Hermitian kernels (matrices).

In the discrete (finite) case, a standard exact sampling algorithm is
introduced in [16] for DPPs with Hermitian kernels. This algorithm uses the
fact that any Hermitian DPP is in fact a mixture of projection DPPs
(elementary DPPs), and also that projection DPPs have a simple exact sampling
algorithm. Other sampling algorithms were also studied, for example see [7,
22], but mostly for DPPs with Hermitian kernels.

However recently a new “greedy” type algorithm was introduced [24, 27], based
on the successive computations of conditional probabilities111Such conditional
approaches were already considered earlier, e.g., [5]. through the block LU
decomposition. In each step of this algorithm we determine whether a given
index is included in the sample or not by a Bernoulli trial, which we refer to
as the observation. Depending on the observation at each step, we modify (or
keep) a diagonal entry (the pivot of the LU decomposition), then perform a
single step of the LU decomposition. This can be less efficient than the
standard Hermitian sampler but this algorithm allows one to sample from non-
Hermitian DPPs.

This paper is inspired by this greedy type algorithm. We make the following
three important remarks:

  * •

After each observation we obtain a kernel corresponding to the new DPP of
unobserved points conditioned on the result of observed points.

  * •

Theoretically, we could force specific points to be included in the sample by
just multiplying Bernoulli parameters, instead of random Bernoulli trials. We
still obtain the above (conditional) kernel under such an event.

  * •

The observations could be done in an arbitrary order (of point indices).

Based on these points, we can use the conditional probability approach to
derive several new determinantal expressions of various probability density
functions (PDF) and cumulative distribution functions (CDF) in Section 3.

Not to be understated is the role of algorithms from numerical linear algebra
in this work as an inspiration and algorithmic enhancement. One way or
another, the kernels of the DPP may undergo a matrix factorization; A
potential key to effective algorithms is the recognition of which choice to
use when.

###  1.3. Main tool

From the third remark above, let us imagine observing a specific point $s$
first. Then from the second remark, force an eigenvalue at (an infinitesimal
interval around) $s$. Finally using the first remark, we introduce the
following Proposition which is the key to our results.

######  Proposition 1.1.

Let $K$ be a kernel of an integral operator222The kernel $K:J\times
J\to\mathbb{C}$ is the kernel of some integral operator $\tilde{K}$ on
$L^{2}(J)$ as follows:  $\tilde{K}f(x)=\int K(x,y)f(y)dy.$ However we will
simply denote by $K$ both the kernel and the integral operator since there is
no confusion throughout this work. on $J$ that defines a continuous DPP of the
$n$-point correlation function (1.1) of some random matrix eigenvalues. Fix a
point $s\in J$ and define a derived kernel

(1.4) $K^{(s)}(x,y):=K(x,y)-\frac{K(x,s)K(s,y)}{K(s,s)}.$

Then, the $n$-point correlation function $p^{(s)}(x_{1},\dots,x_{n})$ of the
rest of eigenvalues given that an eigenvalue already exists in an
infinitesimal interval around $s$ is

$p^{(s)}(x_{1},\dots,x_{n})=\det\left(\left[K^{(s)}(x_{i},x_{j})\right]_{i,j=1,\dots,n}\right).$

In other words, the kernel $K^{(s)}$ defines a continuous DPP (1.1) of the
eigenvalues conditioned on the event of an eigenvalue existing at $s$.

The concept of the conditional DPP in Proposition 1.1 could be found (in terms
of $L$-ensemble) and justified as a DPP in Borodin and Rains [5], where it is
proven to be useful when proving the Eynard-Mehta theorem. It has also been
discussed in the context of machine learning [22, 23].

One might notice that the kernel (1.4) is in fact the result of a single step
of the LU decomposition with the pivot $K(s,s)$, as discussed in the second
remark. Note that we could condition on the selection of more than one points
(Proposition 3.1). We emphasize that this kernel is easy-to-use since it is
explicit and does not include any infinite summation or differentiation.

Proposition 1.1 leads to some new expressions on eigenvalue statistics of
random matrices. In Section 3 we provide several eigenvalue statistics in
terms of Fredholm determinants, which are known to be amenable to numerical
computation through the method proposed in [1, 2]. These results include but
are not limited to:

  * •

PDF, CDF of the two extreme eigenvalues (Sections 3.2, 3.3)

  * •

Joint PDF, CDF of the $k$ extreme eigenvalues (Section 3.4)

  * •

PDF, CDF of the first eigenvalue spacing (Section 3.4)

In these results, a random matrix could be chosen to be any random matrices
with determinantal $n$-point correlation function (1.1), such as the GUE, LUE,
JUE, soft-edge scaling, hard-edge scaling, etc.

###  1.4. Preview #1: Joint PDF of the two largest eigenvalues

A good example of our result is the joint PDF $f^{(\lambda_{1},\lambda_{2})}$
of the $k=2$ largest eigenvalues333The choice of $k=2$ is arbitrary and for
illustrative purpose. We could obtain a joint PDF of any $k$ largest
eigenvalues which we discuss in Section 3.4. $\lambda_{1}\geq\lambda_{2}$. Let
us use the soft-edge scaling limit of the GUE as an example. It is expressible
in terms of a Fredholm determinant444Throughout this paper, we again abuse the
notation for simplicity, denoting the restriction
$K\\!\\!\restriction_{L^{2}(J)}$ of the operator $K$ to square integrable
functions $L^{2}(J)$ by $K\\!\\!\restriction_{J}$. using Proposition 1.1,

(1.5)
$f^{(\lambda_{1},\lambda_{2})}(x_{1},x_{2})=\det\left(\begin{bmatrix}K(x_{1},x_{1})&K(x_{1},x_{2})\\\
K(x_{2},x_{1})&K(x_{2},x_{2})\end{bmatrix}\right)\cdot\det(I-K^{(x_{1},x_{2})}\\!\\!\restriction_{(x_{2},\infty)}),$

for $x_{1}>x_{2}$ and vanishes otherwise, where $K=K_{\operatorname{Ai}}$ is
the Airy kernel (3.2) and the kernel $K^{(x_{1},x_{2})}$ is defined in terms
of $K$,

$K^{(x_{1},x_{2})}(x,y):=K(x,y)-\begin{bmatrix}K(x,x_{1})\\\
K(x,x_{2})\end{bmatrix}^{T}\\!\\!\begin{bmatrix}K(x_{1},x_{1})&K(x_{1},x_{2})\\\
K(x_{2},x_{1})&K(x_{2},x_{2})\end{bmatrix}^{-1}\\!\begin{bmatrix}K(x_{1},y)\\\
K(x_{2},y)\end{bmatrix}.$

Using the above formula (1.5) we were able to compute the correlation
coefficient of the two largest eigenvalues at the soft-edge scaling limit

$\rho(\lambda_{1},\lambda_{2})=0.505\,647\,231\,59...,$

up to 11 digits in less than 2 minutes. This correlation coefficient has a
previously reported computing time of 16 hours for 11 digits in 2010 [1]. The
formula could also be generalized to the $k$ largest or, similarly, smallest
eigenvalues of any random matrices when the $n$-point correlation function of
eigenvalues is given in determinantal manner (1.1). See Section 3.4 for
details.

###  1.5. Preview #2: Determinantal expression for the PDF of the Tracy–Widom
distribution

The famous Tracy–Widom distribution PDF is often plotted. Interestingly, as
far as we know555Except for a recent approach suggested in [3]. See Section
3.2 for details., the last step of the computation involves taking the
derivative of the CDF. In this preview we give a direct determinantal
expression of the Tracy–Widom PDF that give the plot in the bottom left part
of Figure 1. The CDF expression that is typically used is
$F_{2}(s)=\det(I-K_{\operatorname{Ai}}\\!\\!\restriction_{(s,\infty)})$. Our
idea is that we first fix a level at $s$ and then compute the probability that
nothing lies above $s$ with the conditional DPP kernel
$K_{\operatorname{Ai}}^{(s)}$.

######  Proposition 1.2 (PDF of the Tracy–Widom distribution).

The probability density function $f_{2}$ of the largest eigenvalue at the
soft-edge is

(1.6)
$\frac{d}{ds}F_{2}(s)=f_{2}(s)=K_{\operatorname{Ai}}(s,s)\cdot\det\left(I-K_{\operatorname{Ai}}^{(s)}\\!\\!\restriction_{(s,\infty)}\right),$

where $F_{2}$ is the Tracy–Widom distribution (CDF), $K_{\operatorname{Ai}}$
is the Airy kernel and

$K_{\operatorname{Ai}}^{(s)}(x,y)=K_{\operatorname{Ai}}(x,y)-\frac{(\operatorname{Ai}(x)\operatorname{Ai}^{\prime}(s)-\operatorname{Ai}(s)\operatorname{Ai}^{\prime}(x))(\operatorname{Ai}(s)\operatorname{Ai}^{\prime}(y)-\operatorname{Ai}(y)\operatorname{Ai}^{\prime}(s))}{(x-s)(y-s)(s\operatorname{Ai}(s)^{2}-\operatorname{Ai}^{\prime}(s)^{2})}$

is the derived kernel of the conditional DPP as proposed in (1.4).

This again is not limited to the soft-edge scaling limit, but also applicable
to other random matrices such as the finite GUE, LUE, hard-edge and more. See
Section 3.2 for further details. We also provide numerical experiments.

###  1.6. Outline of the paper

In Section 2, we review the theory of discrete and continuous DPPs and their
sampling algorithms. We propose a hybrid sampling method, Algorithm 4, for
DPPs with non-Hermitian projection kernels, for example the DPP of the Aztec
diamond [20]. We then demonstrate an efficient sampling of the DR paths (see
Section 2.3) without sampling the whole Aztec diamond.

Random matrix applications are discussed in Section 3. In Section 3.1 we
review some basic random matrix eigenvalues and the conditional DPP approach.
Then we derive several new determinantal representations which are efficiently
implemented for numerical computations in later sections. In Section 3.2 we
obtain a Fredholm determinant expression for the PDF of extreme eigenvalues,
such as the Tracy–Widom distribution. In Section 3.3 we specialize on the
second largest eigenvalue and provide new formulae for distribution functions.
Section 3.4 discusses the joint PDF of the $k$ largest (extreme) eigenvalues.
Applications of these joint PDFs include the first eigenvalue spacing,
correlation coefficient of the two largest eigenvalues, and many more. Finally
in Section 3.5 we demonstrate the sampling of Airy processes from the DPP.
Throughout Section 3 we provide extensive numerical experiments, and all codes
could be found online.

##  2\. Determinantal point processes

###  2.1. Discrete and Projection DPPs

Discrete DPPs have a fairly straightforward definition as we saw from (1.3).
In particular if we have a finite sized ground set $G$, the kernel is a
(finite) $|G|\times|G|$ matrix $K$, often called the marginal kernel. For a
given DPP $\mathcal{J}$ the following identities are important:

(2.1) $\displaystyle\operatorname{tr}(K)=\mathbb{E}(|\mathcal{J}|),$ (2.2)
$\displaystyle\det(I-K)=\mathbb{P}(\text{$|\mathcal{J}|=0$}).$

It is known that if the kernel matrix $K$ is a projection matrix, a DPP has
the following special property: a DPP with a rank $r$ projection marginal
kernel draws a sample with the size exactly $r$, i.e. $|J|=r$ (almost surely
when continuous).

This property is a cornerstone of the sampling algorithm introduced in [16].
Imagine an algorithm that draws samples from a given DPP, where sample points
are selected in an unsorted (uniformly permuted) order. The probability
$P_{i}$ that a given index $i\in G$ is picked ‘first’ is the following.

$P_{i}=\mathbb{P}(\text{$\mathcal{J}=\\{i\\}$})+\sum_{\begin{subarray}{c}i\in
S\\\
|S|=2\end{subarray}}\frac{1}{2}\mathbb{P}(\text{$\mathcal{J}=S$})+\sum_{\begin{subarray}{c}i\in
S\\\
|S|=3\end{subarray}}\frac{1}{3}\mathbb{P}(\text{$\mathcal{J}=S$})+\cdots=\sum_{i\in
S}\frac{1}{|S|}\mathbb{P}(\text{$\mathcal{J}=S$}).$

Thinking the other way around, if we know $\\{P_{i}\\}$, we could sample the
first point (without worrying about additional sample points) according to the
discrete random variable $X$ defined by $\mathbb{P}(\text{$X=i$})=P_{i}$.
However the probabilities $\\{P_{i}\\}$ in general do not have a simpler
expression in terms of the entries of the marginal kernel.

Nonetheless, for projection DPPs, (normalized) diagonal entries of $K$ equals
the probabilities $\\{P_{i}\\}$ as $P_{i}=K_{ii}/r$ is deduced from (2.1)
together with $\operatorname{tr}(K)=r$. Thus, sampling from a projection DPP
begins by drawing a single index point from a categorical random variable with
normalized diagonal entries as its distribution. After drawing a first point,
one could modify the kernel so we could sample points iteratively as we
describe in the following paragraphs.

Algorithm 1 OrthoProjDPP: Sample from an orthogonal projection DPP

function OrthoProjDPP(Y) $\,\,\,\%\,\,Y\in\mathbb{R}^{n\times r}$ is
orthogonal and $K=YY^{T}$

sample $\leftarrow$ empty vector

for $i=1:r$ do

Draw $j$ from $\\{1,\dots,n\\}$ with
$\mathbb{P}(\text{$X=j$})=\verb|norm|(Y[j,:])^{2}$

Add $j$ to sample

$Q\leftarrow\verb|Householder|(Y[j,:])$

$Y\leftarrow(YQ)[:,2:\text{end}]$

end for

return sample

end function

Let us for a moment restrict our projection matrix $K\in\mathbb{R}^{n\times
n}$ to be an orthogonal projection matrix, so that we have $K=YY^{T}$ for some
orthogonal (unitary, if complex) matrix $Y\in\mathbb{R}^{n\times r}$. The
probability $P_{i}$ above is then equivalent to the squared norm of the
$i^{\text{th}}$ row of $Y$, $\sum_{j=1}^{r}Y_{ij}^{2}$, divided by $r$. We
multiply a Householder reflector $Q$ [9] of the $i^{\text{th}}$ row of $Y$ on
the right side of $Y$, so that $YQ$ has the $i^{\text{th}}$ row
$(P_{i},0,\dots,0)$. If we let $\tilde{Y}=YQ$ we have

$\mathbb{P}(\text{$j\in\mathcal{J}|i\in\mathcal{J}$})=\sum_{k=2}^{r}\tilde{Y}_{jk}^{2},$

which means that the matrix $Z\in\mathbb{R}^{n\times(r-1)}$ obtained by
deleting the first column of $\tilde{Y}$ (which is again orthogonal) plays the
same role as $Y$ when drawing the first index point $i$. In other words,
$ZZ^{T}$ is the rank $r-1$ marginal kernel of the DPP conditioned on the first
sample index point $i$. Recursively applying this procedure $r$ times, we
obtain Algorithm 1 for orthogonal projection DPPs.

The algorithm introduced by Hough et al. in [16] samples from a Hermitian DPP
using the fact that it is a mixture of projection DPPs via the
eigendecomposition of $K$. Algorithm 2 outlines this sampling algorithm.

Algorithm 2 HermDPP: Sample from a Hermitian DPP

function HermDPP($X,\Lambda$) $\%\,\,$Eigendecomposition $K=X\Lambda X^{T}$,
$K\in\mathbb{R}^{n\times n}$

mask $\leftarrow$ empty vector

for $i=1:n$ do

if Bernoulli($\Lambda[i]$)==1 then

Add $i$ to mask

end if

end for

$Y\leftarrow X[:,\text{mask}]$

return `OrthoProjDPP(`$YY^{T}$`)`

end function

###  2.2. Sampling with conditional probabilites

Algorithm 2 and the algorithms followed thereafter selects a single sample
point of $\mathcal{J}$ at each time step. Thus, the set of points that are
‘not selected’ is determined at the final step of the algorithm. On the other
hand, some recent work [24, 27] uses a different approach based on conditional
probabilities. At each step, rather than sampling an index point, we ‘observe’
a single index point and determine whether it is going to be drawn or not.

A central idea comes from the following block LU decomposition where the rows
and columns are partitioned by $(m,n-m)$,

$K=\begin{bmatrix}I_{m}&0\\\
K_{21}K_{11}^{-1}&K_{22}-K_{21}K_{11}^{-1}K_{12}\end{bmatrix}\begin{bmatrix}K_{11}&K_{12}\\\
0&I_{n-m}\end{bmatrix}.$

This is equivalent to the result of $m$ steps of the (unpivoted) LU
decomposition. We have the following conditional probability for any subset
$S$ of $\\{m+1,\dots,n\\}$,

$\det\left(\left[(K_{22}-K_{21}K_{11}^{-1}K_{12})_{i,j}\right]_{i,j\in
S}\right)=\mathbb{P}(\text{$S\subset\mathcal{J}$ $|$
$\\{1,\dots,m\\}\subset\mathcal{J}$}).$

In other words, the matrix $K_{22}-K_{21}K_{11}^{-1}K_{12}$ (i.e., the Schur
complement) serves as a new kernel for the DPP conditioned on all of
$1,\dots,m$ already being drawn. Recall from the remark in the introduction
that indices $1,\dots,m$ could be in fact any $m$ indices, since the order of
observation could be arbitrarily chosen by some row/column permutation.
Furthermore, a DPP with the condition that some index points are not selected
could also be analyzed in a similar manner. The following Proposition
summarizes this idea.

######  Proposition 2.1 ([27]).

Let $K$ be the kernel of the DPP $\mathcal{J}$. Given disjoint subsets $A,B$
of the ground set $\\{1,\dots,n\\}$, we have the following probabilities.

$\displaystyle\mathbb{P}(\text{$B\subset\mathcal{J}$ $|$
$A\subset\mathcal{J}$})$
$\displaystyle=\det(K_{B,B}-K_{B,A}K_{A,A}^{-1}K_{A,B}),$
$\displaystyle\mathbb{P}(\text{$B\subset\mathcal{J}$ $|$
$a\notin\mathcal{J},\forall a\in A$})$
$\displaystyle=\det(K_{B,B}-K_{B,A}(K_{A,A}-I)^{-1}K_{A,B}),$

where $K_{X,Y}$ is the submatrix of $K$ with row indices $X$ and column
indices $Y$.

Performing Bernoulli trials on pivots and using Proposition 2.1 we have the
following Algorithm 3 from [24, 27].

Algorithm 3 genDPP: Sample from a general DPP

function genDPP($K$) $\%\,\,K\in\mathbb{R}^{n\times n}$

sample $\leftarrow$ empty vector

for $i=1:n$ do

if Bernoulli($K[i,i]$)==1 then

Add $i$ to sample

else

$K[i,i]\leftarrow K[i,i]-1$

end if

$K[i+1:n,i+1:n]\,\,-=\,\,K[i+1:n,i]*K[i,i+1:n]/K[i,i]$

end for

return sample

end function

Although this “greedy” type approach may be inefficient, there is one
significant advantage: Algorithm 3 enables sampling of a non-Hermitian DPP.
For example the discretized DPP of Dyson Brownian motion, which we will
discuss in Section 3.5, is a non-Hermitian DPP. Another example is the Aztec
diamond domino tiling [17, 18, 20] which is used in [27] as an example of a
non-Hermitian DPP.

Algorithm 4 nonOrthoProjDPP: Sample from a non-Hermitian projection DPP

function nonOrthoProjDPP($K$) $\%\,\,K\in\mathbb{R}^{n\times n}$ non-Hermitian
rank $r$ projection

sample $\leftarrow$ empty vector

for $i=1:r$ do

Draw $j$ from $\\{1,\dots,n\\}$ with $\mathbb{P}(\text{$X=j$})=K[j,j]/(r-i+1)$

Add $j$ to sample

$K\,\,-=\,\,K[:,j]*K[j,:]/K[j,j]$

end for

return sample

end function

However one might notice that the DPP kernel of the Aztec diamond obtained
from Kenyon’s formula using the inverse Kasteleyn matrix [6, 21] is non-
Hermitian but still a projection matrix. (A clue would be that any sample
always has the fixed size $n(n+1)=\text{number of dominos}$.) When we have a
non-Hermitian projection DPP with a small rank, Algorithm 3 could be improved
by combining with Algorithm 1. At each step we draw indices as we did in
Algorithm 1 from diagonal entries and then we modify the kernel as in
Algorithm 3. This reduces the number of steps in Algorithm 3 from the size of
$K$ to its rank. We briefly sketch this hybrid approach in Algorithm 4.

We additionally note that unfortunately, the Aztec diamond DPP has its rank
and size of the same order, which only obtains a small or even no improvement
in practice. Nonetheless, if a kernel $K$ is non-Hermitian projection with
$\text{rank}(K)\ll n$, theoretically Algorithm 4 should outperform Algorithm
3. Table 1 is the summary of appropriate choices of exact DPP samplers,
depending on the kernel $K$.

Alg # | Hermitian? | Projection? | Example
---|---|---|---
1 | $\color[rgb]{0,0.45,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.45,0}\checkmark$ | $\color[rgb]{0,0.45,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.45,0}\checkmark$ | Hermite kernel (finite GUE)
2 | $\color[rgb]{0,0.45,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.45,0}\checkmark$ | x | Airy kernel (soft-edge truncated)
3 | x | x | Airy process (Section 3.5)
4 | x | $\color[rgb]{0,0.45,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.45,0}\checkmark$ | Aztec diamond
Table 1. Choice of an exact sampling algorithm depending on properties of the
kernel of a DPP.

###  2.3. Application to Aztec diamonds: efficiently sampling a DR path

Not only can Algorithm 3 sample non-Hermitian DPPs but also it enables partial
sampling. One application that benefits from the partial sampling of Algorithm
3 is the sampling of DR paths [30] of the Aztec diamond. The north polar
region (NPR) boundary process is of a great interest due to its connection to
corner growth, and eventually to the Airy process [20]. DR paths of the Aztec
diamond are defined as follows: We label the vertical dominos east-going (E)
and west-going (W), according to their checkerboard patterns. An easy way to
remember is that the vertical domino that fits the westmost corner is a
W-domino. Similarly we label horizontal dominos N-dominos and S-dominos. In
the left of Figure 2, red, blue, green, yellow dominos are W, E, N, S-dominos,
respectively. We draw horizontal lines in the middle of S-dominos, $\pm 45$
degree lines passing through the center of W, E-dominos, respectively. The
right of Figure 2 has such lines drawn in red. It is known that these lines
form $n$ continuous paths, which are called the DR paths.

Figure 2. An Aztec diamond domino tiling with $n=10$ sampled from a DPP (left)
and corresponding DR paths in red (right)

An interesting observation is that we do not need the whole Aztec domino
configuration to get the top DR path. Using Algorithm 3 partially we could
efficiently sample the top DR path (and similarly other specific DR paths) by
only observing the possible dominos along the path. For example let us start
by observing the westmost location. If we sample a W-domino there we then
‘observe’ the three next possible (W, E, S) dominos that share the upper half
of the right side of the westmost domino, since they are the three possible
connections to the current path. This is possible due to the last point of the
remark in the introduction, that we could observe in any desired order. We
recursively do this until we reach the east boundary. This partial sampling
along the top DR path reduces the complexity of DPP sampling by $O(n)$. Figure
3 illustrates this point by comparing the two DPP sampling times, where
sampling only the top DR path is about 40 times faster.

DPP Sample time: 52.44 sec DPP Sample time: 1.28 sec

Figure 3. An Aztec diamond domino tiling with $n=30$ sampled from a DPP (left)
and a top DR path sampled without the whole domino configuration (right). The
sampling time for the whole domino tiling is 52.44 seconds and the sampling
time for a top DR path is 1.28 second. Moreover, with 50 seconds one could
sample a top DR path with $n=50$. Nevertheless, we have not found this
algorithm to be yet competitive with usual Aztec domino sampling algorithms,
e.g., the shuffling algorithm [10].

###  2.4. Continuous DPPs

Many concepts in discrete DPPs extend to continuous DPPs. The ground set $G$
now becomes a continuous interval (or any set) and the marginal kernel matrix
becomes the kernel function $K:G\times G\to\mathbb{C}$. A continuous DPP
defines the following $n$-point correlation function, which is the continuous
analogue of $\mathbb{P}(\text{$S\subset\mathcal{J}$})$ in (1.3),

$p(x_{1},\dots,x_{n})=\det\left([K(x_{i},x_{j})]_{i,j=1,\dots,n}\right).$

One way to describe the $n$-point correlation function $p(x_{1},\dots,x_{n})$
is the following:

(2.3) $\lim_{\delta x\to 0}\frac{1}{(\delta x)^{n}}{\mathbb{P}(\text{$n$
points at length $\delta x$ intervals around $x_{1},\dots,x_{n}$})}.$

Sampling algorithms for a continuous DPP could be generalized from the
discrete case. For a continuous projection DPP with bounded trace, i.e.,

$\int K(x,y)K(y,z)dy=K(x,z)\hskip 14.22636pt\text{and}\hskip 14.22636pt\int
K(x,x)dx<\infty,$

Algorithm 4 could be generalized easily. The point selection at each timestep
is now a univariate random variable with its PDF proportional to the diagonal
$K(x,x)$. More details could be found in [15].

Another method for sampling from a continuous DPP is by discretizing a
continuous DPP to a discrete DPP. Let us use the Hermite kernel
$K_{\text{Herm}}^{(N)}$, (1.2) as an example. To create a finite matrix we
truncate the ground set. For the Hermite kernel with $N=5$, it is known that
the largest eigenvalue lies around $\sqrt{2N}\approx 3$, and the probability
that an eigenvalue lies outside a wide interval of length $L$, for example
$(-10,10)$, is already much lower than the machine epsilon of the double
precision. Then with $M\approx L/\delta x$ length $\delta x$ intervals in the
truncated region, create an $M\times M$ matrix $K$ such that
$K_{ij}=K(x_{i},x_{j})\delta x$, where $x_{j}$ is the midpoint of the
$j^{\text{th}}$ interval.666Of course, as $\delta x\to 0$ this does not have
to be the midpoint. From (2.3) the determinants of principal submatrices
approach probabilities of sample points lying around corresponding intervals,
as $\delta x\to 0$. On the other hand if one tries to compute integrals such
as $\operatorname{tr}K=\int K(x,x)dx$ or the Fredholm determinant $\det(I-K)$,
one may use weights and points from quadrature rules, see [2] for details. In
Section 3.5 we discretize the extended Airy kernel to sample Airy processes
using Algorithm 3.

##  3\. The conditional DPP method applied to random matrix theory

Let $K$ be a kernel for an integral operator such that the $n$-point
correlation function $p$ of some random matrix (or its scaling limit)
eigenvalues is given as an $n\times n$ determinant (1.1). Examples of such
random matrices are the GUE, LUE, soft-edge, hard-edge, etc. As mentioned
above, $p(x_{1},\dots,x_{n})$ is the continuous analogue of the set-level
complementary cumulative distribution function (CCDF),
$\mathbb{P}(\text{$S\subset\mathcal{J}$})$ in (1.3).

The probability that no eigenvalue exists in a given interval $J$ plays an
important role. In the continuous case, it is given as the following Fredholm
determinant

(3.1) $\mathbb{P}(\text{No eigenvalue in
$J$})=\det(I-K\\!\\!\restriction_{J}),$

which is a continuous generalization of (2.2).

###  3.1. Kernel of the DPP from conditional probabilities

Let us first prove Proposition 1.1, which is the continuous analogue of
Proposition 2.1.

###### Proof of Proposition 1.1.

One step of the LU decomposition (in reverse order) of the $(n+1)$-point
correlation function is

$[K(x_{i},x_{j})]_{i,j=1,\dots,n+1}=\begin{bmatrix}\left[K^{(x_{n+1})}(x_{i},x_{j})\right]_{i,j=1,\dots,n}&v\\\
0&1\end{bmatrix}\begin{bmatrix}I_{n}&0\\\ w&K(x_{n+1},x_{n+1})\end{bmatrix},$

where $v=[K(x_{i},x_{n+1})]_{i=1,\dots,n}/K(x_{n+1},x_{n+1})$ and
$w^{T}=[K(x_{n+1},x_{i})]_{i=1,\dots,n}$. From the multiplicativity of
determinants we get

$\det\left(\left[K(x_{i},x_{j})\right]_{i,j=1,\dots,n+1}\right)=\det\left(\left[K^{(x_{n+1})}(x_{i},x_{j})\right]_{i,j=1,\dots,n}\right)\cdot
K(x_{n+1},x_{n+1}).$

Then, since
$p^{(x_{n+1})}(x_{1},\dots,x_{n})=\det([K^{(x_{n+1})}(x_{i},x_{j})]_{i,j=1,\dots,n})$
and using (2.3) for the left hand side and $K(x_{n+1},x_{n+1})$ we have

$\displaystyle p^{(x_{n+1})}$ $\displaystyle(x_{1},\dots,x_{n})$
$\displaystyle=\lim_{\delta x\to 0}\frac{1}{(\delta
x)^{n}}\frac{\mathbb{P}(\text{Eigenvalues at $\delta x$ intervals around
$x_{1},\dots,x_{n+1}$})}{\mathbb{P}(\text{Eigenvalue at $\delta x$ interval
around $x_{n+1}$})},$

which is the desired $n$-point correlation function conditioned on the event
that an eigenvalue existing around an infinitesimal interval containing
$x_{n+1}$. ∎

Note that we could also generate a kernel for an $n$-point correlation
function conditioned on the event that a fixed location does not contain an
eigenvalue, by replacing the denominator of the kernel (1.4) with $K(s,s)-1$.
However, we note that the condition that an eigenvalue does not exist at a
specific location is less versatile for applications than the condition that
an eigenvalue exists at a specific location, in the continuous setting.

We generalize Proposition 1.1 to a conditional DPP kernel with any number
($m$) of forced eigenvalues.

######  Proposition 3.1.

With the same assumptions as in Proposition 1.1, the $n$-point correlation
function $p^{(s_{1},\dots,s_{m})}$ of the eigenvalues conditioned on the event
that $m$ eigenvalues already exist around infinitesimal intervals around
$s_{1},\dots,s_{m}$ is

$p^{(s_{1},\dots,s_{m})}(x_{1},\dots,x_{n})=\det\left(\left[K^{(s_{1},\dots,s_{m})}(x_{i},x_{j})\right]_{i,j=1,\dots,n}\right),$

where the kernel $K^{(s_{1},\dots,s_{m})}$ is given as

$\displaystyle K$ ${}^{(s_{1},\dots,s_{m})}(x,y)$
$\displaystyle=K(x,y)-\begin{bmatrix}K(x,s_{1})\\\ \vdots\\\
K(x,s_{m})\end{bmatrix}^{T}\\!\\!\begin{bmatrix}K(s_{1},s_{1})&\cdots\\!\\!&K(s_{1},s_{m})\\\
\vdots&\ddots\\!\\!&\vdots\\\
K(s_{m},s_{1})&\cdots\\!\\!&K(s_{m},s_{m})\end{bmatrix}^{-1}\\!\\!\begin{bmatrix}K(s_{1},y)\\\
\vdots\\\ K(s_{m},y)\end{bmatrix}.$

###### Proof.

The proof is analogous to the proof of Proposition 1.1 replacing the LU
decomposition with the block LU decomposition of an $(n+m)\times(n+m)$ matrix,
with row/column partitions $(n,m)$. ∎

###  3.2. The PDF of extreme (largest) eigenvalues

One can derive a new expression for the PDF of the largest eigenvalue of a
random matrix. Let us first review some basic eigenvalue statistics and
conventions. For the $N\times N$ GUE we have

$E_{2}^{(N)}(0;J):=\mathbb{P}(\text{No eigenvalue of the $N\times N$ GUE is in
$J$})=\det(I-K_{\text{Herm}}^{(N)}\\!\\!\restriction_{J}),$

where the Hermite kernel is defined in (1.2). In particular we denote
$E_{2}^{(N)}(0;(0,s))$ simply by $E^{(N)}(0;s)$.

As $N\to\infty$ some scaling limits are defined. With the sine kernel
$K_{\sin}(x,y)=\frac{\sin\pi(x-y)}{\pi(x-y)}$ one obtains the bulk scaling
limit777In fact, appropriate scaling of $\beta=2$ Laguerre and Jacobi
ensembles in the bulk also yield the same bulk limit [26]. defined with mean
spacing 1,

$E(0;s):=\lim_{N\to\infty}E^{(N)}(0;\frac{\pi}{\sqrt{2n}}s)=\det\left(I-K_{\sin}\\!\\!\restriction_{(0,s)}\right).$

Also with the Airy kernel

(3.2)
$K_{\operatorname{Ai}}(x,y)=\frac{\operatorname{Ai}(x)\operatorname{Ai}^{\prime}(y)-\operatorname{Ai}^{\prime}(x)\operatorname{Ai}(y)}{x-y},$

we have the following Fredholm determinant representation of the soft-edge
scaling limit of the GUE888As in the bulk, this could also be the largest
eigenvalue of the LUE in the appropriate soft-edge scaling limit.[12] (the
Tracy–Widom distribution $F_{2}$)

$F_{2}(s):=\\!\lim_{N\to\infty}\\!E_{2}^{(N)}\\!\\!\left(0;(\sqrt{2n}+\frac{s}{\sqrt{2n^{1/3}}},\infty)\right)\\!\\!=\mathbb{P}(\text{$\lambda_{\text{max}}\leq
s$})=\det\left(I-K_{\operatorname{Ai}}\\!\\!\restriction_{(s,\infty)}\right).$

Depending on the position and restrictions imposed on eigenvalues, $\det(I-K)$
turns into several different distribution functions. At the edge (either hard
or soft) $\det(I-K)$ becomes the CDF. In particular, at the $+\infty$ side of
the soft-edge, $\det(I-K_{\operatorname{Ai}}\\!\\!\restriction_{(s,\infty)})$
equals $\mathbb{P}(\text{$\lambda_{\text{max}}\leq s$})$, the CDF of the
largest eigenvalue. Similarly, at the LUE hard-edge999Also Jacobi ensemble
with $\beta=2$ has the same hard-edge scaling limit.[4]
$\det(I-K_{\operatorname{Bess},\alpha}\\!\\!\restriction_{(0,s)})$ equals
$\mathbb{P}(\text{$\lambda_{\text{min}}\geq s$})$, the CCDF.

On the other hand in the bulk, $\det(I-K)$ is not a CDF. Rather, its
derivative is a CDF, i.e., $E(0;s)$ in the bulk has the following derivatives
[25, Chapter 6.1.2],

(3.3) $\tilde{F}(0;s):=-\frac{d}{ds}E(0;s),\hskip
28.45274ptp(0;s):=-\frac{d}{ds}\tilde{F}(0;s).$

The first derivative $\tilde{F}(0;s)$ is the probability (in the bulk) that,
for a randomly chosen level around zero, the distance $s$ to the right
contains no eigenvalue. If we let a random variable $D$ to be the distance
(again, starting from any randomly chosen level) to the right until the next
level, $\tilde{F}(0;s)=\mathbb{P}(D>s)$ is a CCDF. It follows that $p(0;s)$ is
the PDF of $D$.

Some of these probabilities can alternatively be described by conditional
probabilities. For example, the PDF of the smallest eigenvalue at the hard-
edge is the product of (1) the $1$-point correlation function at $s$
($=K_{\operatorname{Bess},\alpha}(s,s)$) and (2) the probability that no
eigenvalue lies in $(0,s)$, conditioned on an eigenvalue existing around $s$.
The latter could be obtained from Proposition 1.1 and becomes a new expression
for the PDF of the smallest eigenvalue at the hard-edge,

(3.4)
$f_{\text{hard},\alpha}(s)=\underbrace{K_{\operatorname{Bess},\alpha}(s,s)}_{\text{(1)
Level at $s$}}\hskip
5.69046pt\cdot\underbrace{\det\left(I-K_{\operatorname{Bess},\alpha}^{(s)}\\!\\!\restriction_{(0,s)}\right)}_{\text{(2)
No levels on $(0,s)$, given level at $s$}}.$

In Proposition 1.2 we have already introduced the PDF of the Tracy–Widom
distribution with the same idea. Generalizing these, we get the following
Corollary.

######  Corollary 3.1.1.

Given a kernel $K$ as in Proposition 1.1 that it is a kernel of a continuous
DPP of some random matrix eigenvalues. Let

$f(a,b)=\det(I-K\\!\\!\restriction_{(a,b)}),$

the probability that no eigenvalue lies in $J=(a,b)$. Then, the following
holds:

(3.5) $\displaystyle\frac{d}{da}f(a,b)$
$\displaystyle=K(a,a)\det(I-K^{(a)}\\!\\!\restriction_{J})=\mathbb{P}(\text{eigenvalue
at $a$ and none in $J$}),$ (3.6) $\displaystyle-\frac{d}{db}f(a,b)$
$\displaystyle=K(b,b)\det(I-K^{(b)}\\!\\!\restriction_{J})=\mathbb{P}(\text{eigenvalue
at $b$ and none in $J$}),$

with the kernels $K^{(a)}$ and $K^{(b)}$ defined as in (1.4).

###### Proof.

We prove (3.5) and the proof for (3.6) is similar. Let $a^{\prime}=a+\delta
a$.

$\displaystyle\frac{1}{\delta a}($ $\displaystyle
f(a^{\prime},b)-f(a,b))=\frac{1}{\delta a}\left(\mathbb{P}(\text{Nothing in
$(a^{\prime},b)$})-\mathbb{P}(\text{Nothing in $(a,b)$})\right)$
$\displaystyle=\delta a^{-1}\mathbb{P}(\text{Eigenvalue at $(a,a^{\prime})$,
none in $(a^{\prime},b)$})$ $\displaystyle=\delta a^{-1}\mathbb{P}(\text{None
in $(a^{\prime},b)\,|$ Eigenvalue at
$(a,a^{\prime})$})\cdot\mathbb{P}(\text{Eigenvalue at $(a,a^{\prime})$})$
$\displaystyle=\delta
a^{-1}\det(I-K^{(a)}\\!\\!\restriction_{(a^{\prime},b)})\cdot\mathbb{P}(\text{Eigenvalue
at $(a,a^{\prime})$})$

As we let $\delta a\to 0$, $\delta a^{-1}\mathbb{P}(\text{Eigenvalue at
$(a,a^{\prime})$})$ goes to $K(a,a)$. ∎

Let us also apply Corollary 3.1.1 to distribution functions (3.3) in the bulk.
For $\tilde{F}(0;s)$, since $K_{\sin}(0,0)=1$ we have

(3.7) $\tilde{F}(0;s)=-\det(I-K_{\sin}^{(0)}\\!\\!\restriction_{(0,s)}),$

where (1.4) defines
$K_{\sin}^{(0)}=\frac{\sin\pi(x-y)}{\pi(x-y)}-\frac{\sin\pi x\sin\pi
y}{\pi^{2}xy}$. Applying once more we get

(3.8)
$p(0;s)=K_{\sin}^{(0)}(s,s)\cdot\det(I-K_{\sin}^{(0)*}\\!\\!\restriction_{(0,s)}),$

where the kernel
$K_{\sin}^{(0)*}(x,y):=K_{\sin}^{(0)}(x,y)-\frac{K_{\sin}^{(0)}(x,s)K_{\sin}^{(0)}(s,y)}{K_{\sin}^{(0)}(s,s)}$
is a twice derived sine kernel. Figure 4 is the plot of $\tilde{F}(0;s)$,
$p(0;s)$ using (3.7), (3.8), respectively.

$\tilde{F}(0;s)$ | $p(0;s)$
---|---
|
Figure 4. Plots of $\tilde{F}(0;s)$ (left) and $p(0;s)$ (right) defined in
(3.3). Numerical computation was done by evaluating Fredholm determinants [2]
in equations (3.7), (3.8).

Moreover, Table 2 is the first four moments of $D$ obtained from the
computation of $p(0;D=s)$ values. See code F0p0 for the implementation
resulting in Figure 4 and Table 2.

Mean | Variance | Skewness | Excess Kurtosis
---|---|---|---
1.0 | 0.179 993 877 691 8 | 0.497 063 620 491 8 | 0.126 699 848 039 9
Table 2. The first four moments of the spacing $D$ near zero in the bulk up to
13 digits. Total computation time is 0.052 second. One could compare this
result to Table 8 of [1].

One benefit of our formulae, such as Proposition 1.2, (3.4), (3.7), (3.8), is
an efficient and accurate numerical computation. We provide numerical
experiment that compares our computation of the Tracy–Widom PDF, (1.6),
$f_{2}(s)=K_{\operatorname{Ai}}(s,s)\det(I-K_{\operatorname{Ai}}^{(s)}\\!\\!\restriction_{(s,\infty)})$
with some other numerical approaches. Numerical results show that our
expressions can be as efficient and potentially more accurate.

One approach for computing Tracy–Widom PDF is the automatic differentiation on
$F_{2}(s)$ combined with the Fredholm determinant computation [2] of
$F_{2}(s)=\det(I-K_{\operatorname{Ai}}\\!\\!\restriction_{(s,\infty)})$.
Another approach has also been suggested recently in [3, Eq (37b)],

$f_{2}(s)=-F_{2}(s)\cdot\operatorname{tr}\big{(}(I-K)^{-1}K^{\prime}\big{)}\\!\\!\restriction_{L^{2}(s,\infty)}=F_{2}(s)\cdot\langle(I-K)^{-1}\operatorname{Ai},\operatorname{Ai}\rangle_{L^{2}(s,\infty)}$

where $K^{\prime}$ is defined as
$K^{\prime}(x,y)=\left(\frac{\partial}{\partial x}+\frac{\partial}{\partial
x}\right)K(x,y)$. See equation (34b) and nearby discussion in [3]. In Table 3,
we compare the accuracy of the computation of $f_{2}(s)$, for different number
of quadrature points and $s$ values, using two different methods: (1) Equation
(37b) of [3] and (2) our approach, (1.6). Comparison for other random matrix
statistics such as in the bulk (sine kernel) or hard-edge (Bessel kernel) for
the values of the derivative of log determinant and comparison against the
automatic differentiation could also be found in the code prime-computation.

| $m=10$ | $m=20$
---|---|---
$s\,\,$ | Eq (37b) of [3] | Eq (1.6) | Eq (37b) of [3] | Eq (1.6)
-4.0 | $3.79\times 10^{-1}$ | $2.76\times 10^{-2}$ | $3.72\times 10^{-7}$ | $5.60\times 10^{-7}$
-3.5 | $2.94\times 10^{-2}$ | $6.13\times 10^{-2}$ | $5.14\times 10^{-8}$ | $8.25\times 10^{-9}$
-3.0 | $2.48\times 10^{-2}$ | $1.68\times 10^{-2}$ | $1.07\times 10^{-8}$ | $4.18\times 10^{-9}$
-2.5 | $8.76\times 10^{-3}$ | $1.89\times 10^{-3}$ | $9.41\times 10^{-11}$ | $1.39\times 10^{-10}$
-2.0 | $1.41\times 10^{-3}$ | $2.55\times 10^{-3}$ | $1.28\times 10^{-10}$ | $2.11\times 10^{-13}$
-1.5 | $2.27\times 10^{-3}$ | $1.70\times 10^{-3}$ | $4.70\times 10^{-11}$ | $3.33\times 10^{-12}$
-1.0 | $5.40\times 10^{-4}$ | $3.87\times 10^{-5}$ | $3.80\times 10^{-11}$ | $1.70\times 10^{-12}$
-0.5 | $5.61\times 10^{-4}$ | $3.06\times 10^{-6}$ | $7.12\times 10^{-12}$ | $1.87\times 10^{-14}$
0.0 | $5.92\times 10^{-4}$ | $1.17\times 10^{-6}$ | $1.15\times 10^{-11}$ | $9.35\times 10^{-14}$
Table 3. Absolute values of relative errors of two approaches for computing
the PDF of the Tracy–Widom distribution $f_{2}(s)$, with $m$ point Gauss-
Legendre quadratures. The expression
$f_{2}(s)=F_{2}(s)\cdot\langle(I-K)^{-1}\operatorname{Ai},\operatorname{Ai}\rangle_{L^{2}(s,\infty)}$
is used from Equation (37b) of [3] and our approach uses (1.6),
$f_{2}(s)=K_{\operatorname{Ai}}(s,s)\cdot\det(I-K_{\operatorname{Ai}}^{(s)}\\!\\!\restriction_{(s,\infty)})$.
Our approach shows slightly better overall accuracy but the errors converge
very quickly to the machine precision in both methods.

###  3.3. The PDF and CDF of the second largest eigenvalue

Some new formulae for the second largest (second extreme) eigenvalue could be
obtained from Proposition 1.1. Let us again take the soft-edge and the Airy
kernel (3.2) as an example. A standard formula for the CDF of the
$k^{\text{th}}$ largest eigenvalue is

(3.9)
$F_{2}(k;s)=\sum_{m=0}^{k-1}\frac{(-1)^{m}}{m!}\frac{d^{m}}{dz^{m}}\det\left(I-zK_{\operatorname{Ai}}\\!\\!\restriction_{(s,\infty)}\right)\big{|}_{z=1}.$

When $k=2$ we could derive a somewhat different formula for the CDF and PDF
that does not involve differentiation or summation using the conditional DPP.
For the CDF $\mathbb{P}(\lambda_{2}<s)$, we need to compute the probability
that there is only one level lying in $(s,\infty)$ when $s$ is given. In the
discrete DPP, the probability of having only one sample point is the trace of
the $L$ kernel, where $L=(I-K)^{-1}K$, divided by $\det(I+L)$, which also
holds similarly in the continuous case. Thus we obtain

$F_{2}(2;s)=\underbrace{\operatorname{tr}\left((I-K_{\operatorname{Ai}})^{-1}K_{\operatorname{Ai}}\right)\\!\\!\restriction_{(s,\infty)}}_{\operatorname{tr}(L)}\cdot\underbrace{\det(I-K_{\operatorname{Ai}}\\!\\!\restriction_{(s,\infty)})}_{\det(I+L)^{-1}}.$

For the PDF of the second largest eigenvalue, we need to fix an eigenvalue at
$s$ and proceed similarly as in the CDF. More precisely, we multiply (1) the
$1$-point correlation function at $s$ and (2) the probability that there is
only a single eigenvalue above $s$, conditioned on (1). That is,

$\displaystyle f_{2}(2;s)$ $\displaystyle=\frac{d}{ds}F_{2}(2;s)$ (3.10)
$\displaystyle=\underbrace{K_{\operatorname{Ai}}(s,s)}_{\text{level at
$s$}}\cdot\underbrace{\operatorname{tr}\left((I-K_{\operatorname{Ai}}^{(s)})^{-1}K_{\operatorname{Ai}}^{(s)}\right)\\!\\!\restriction_{(s,\infty)}\det(I-K_{\operatorname{Ai}}^{(s)}\\!\\!\restriction_{(s,\infty)})}_{\text{Only
one eigenvalue in $(s,\infty)$ given a level at $s$}}.$

Generalizing this, we obtain the following proposition.

######  Proposition 3.2.

For a random matrix that has its eigenvalue $n$-point correlation function
given as in Proposition 1.1 with a kernel $K$, the CDF and PDF of the second
largest eigenvalue $\lambda_{2}$ are

(3.11) $\displaystyle
F^{\lambda_{2}}(s)=\operatorname{tr}\left((I-K)^{-1}K\right)\\!\\!\restriction_{(s,\infty)}\cdot\det(I-K\\!\\!\restriction_{(s,\infty)}),$
(3.12) $\displaystyle
f^{\lambda_{2}}(s)=K(s,s)\cdot\operatorname{tr}\left((I-K^{(s)})^{-1}K^{(s)}\right)\\!\\!\restriction_{(s,\infty)}\cdot\det(I-K^{(s)}\\!\\!\restriction_{(s,\infty)}),$

where the kernel $K^{(s)}$ is given as (1.4).

Proposition 3.2 could be used with the Airy kernel (soft-edge) as well as
finite $N$ kernels like the Hermite kernel. In the hard-edge scaling the
interval $(0,s)$ is used instead of $(s,\infty)$ with the Bessel kernel for
the second smallest eigenvalue. Again, expressions (3.11) and (3.12) yield
highly accurate numerical values using Bornemann’s approach [2]. Table 4 is
the computed first four moments of the two smallest eigenvalues
($\lambda_{1}<\lambda_{2}$) in the hard-edge and soft-edge using PDF formulas
of Section 3.2 and (3.12).

|  | Mean | Variance | Skewness | Kurtosis
---|---|---|---|---|---
Soft | $\lambda_{1}$ | -1.771 087 | 0.813 195 | 0.224 084 | 0.093 448
$\lambda_{2}$ | -3.675 437 | 0.540 545 | 0.125 027 | 0.021 740
Hard | $\lambda_{1}$ | 4.000 000 | 16.000 000 | 2.000 000 | 6.000 000
$\alpha=0$ | $\lambda_{2}$ | 24.362 715 | 140.367 319 | 0.924 147 | 1.225 112
Hard | $\lambda_{1}$ | 10.873 127 | 55.745 139 | 1.320 312 | 2.541 266
$\alpha=1$ | $\lambda_{2}$ | 40.812 203 | 259.898 510 | 0.737 801 | 0.764 990
Hard | $\lambda_{1}$ | 20.362 715 | 124.367 319 | 1.015 815 | 1.461 306
$\alpha=2$ | $\lambda_{2}$ | 60.112 814 | 416.851 440 | 0.622 605 | 0.532 483
Table 4. The first four moments (the last column is the excess kurtosis) of
the two extreme eigenvalues, soft-edge and hard-edge. Computation time is less
than a second for the whole table. See codes cor-coeff-softedge and cor-coeff-
hardedge for the implementation.

###  3.4. The joint distribution of the $k$ largest eigenvalues

In this section we derive an expression for the joint distribution of the $k$
largest (or smallest) eigenvalues in terms of a Fredholm determinant.

In the case of $k=2$, the joint density of the two smallest eigenvalues of the
Laguerre ensemble has been studied in [14]. In addition, the joint
distribution of the two smallest eigenvalues at the hard-edge is obtained in
terms of the solution of a differential equation that resembles Jimbo-Miwa-
Okamoto $\sigma$-form of the Painlevé III. Analogously [32], studies the joint
density of the two largest eigenvalues at the soft-edge, obtaining an
expression in terms of Painlevé II transcendents and isomonodromic components
using hard-to-soft edge transition [4].

With Proposition 3.1, we derive an expression for the joint PDF of the $k$
largest eigenvalues in terms of the Fredholm determinant. Let us take for
example $k=2$ and consider the $N\times N$ GUE. For simplicity let $K$ be the
Hermite kernel (1.2) and $p$ the $N\times N$ GUE eigenvalue $n$-point
correlation function. The $n$-point correlation function after forcing (and
conditioning on) two eigenvalues on $x_{1},x_{2}$ is

$p^{(x_{1},x_{2})}(y_{1},\dots,y_{n})=\det\left(\left[K^{(x_{1},x_{2})}(y_{i},y_{j})\right]_{i,j=1,\dots,n}\right).$

where the kernel $K^{(x_{1},x_{2})}$ is given as

$K^{(x_{1},x_{2})}(x,y):=K(x,y)-\begin{bmatrix}K(x,x_{1})\\\
K(x,x_{2})\end{bmatrix}^{T}\begin{bmatrix}K(x_{1},x_{1})&K(x_{1},x_{2})\\\
K(x_{2},x_{1})&K(x_{2},x_{2})\end{bmatrix}^{-1}\begin{bmatrix}K(x_{1},y)\\\
K(x_{2},y)\end{bmatrix}.$

Then, the joint PDF $f^{(\lambda_{1},\lambda_{2})}$ of the two largest
eigenvalues $\lambda_{1}\geq\lambda_{2}$ is obtained as the following
conditional probability argument assuming $x_{1}>x_{2}$,

$\displaystyle f^{(\lambda_{1},\lambda_{2})}(x_{1},x_{2})$
$\displaystyle=p(x_{1},x_{2})\cdot\mathbb{P}(\text{No other eigenvalues in
$(x_{2},\infty)\,\,|\,\,\lambda_{1}=x_{1},\lambda_{2}=x_{2}$})$
$\displaystyle=\det\left(\begin{bmatrix}K(x_{1},x_{1})&K(x_{1},x_{2})\\\
K(x_{2},x_{1})&K(x_{2},x_{2})\end{bmatrix}\right)\cdot\det(I-K^{(x_{1},x_{2})}\\!\\!\restriction_{(x_{2},\infty)}),$

and $f^{(\lambda_{1},\lambda_{2})}(x_{1},x_{2})$ vanishes when $x_{1}\leq
x_{2}$ .

Similarly for the two smallest eigenvalues of the LUE, one could replace the
interval $(x_{2},\infty)$ of the above right hand side with $(0,x_{2})$, where
$0\leq x_{1}\leq x_{2}$ being the two smallest eigenvalues. Again this
approach works for any random matrix levels with determinantal $n$-point
correlation function.

Extending to general $k$’s we obtain an expression for the joint PDF of the
$k$ largest eigenvalues, $f^{(\lambda_{1},\dots,\lambda_{k})}$, as the
following.

######  Proposition 3.3.

For a random matrix whose eigenvalue $n$-point correlation function is given
in determinantal representation (2.3) with a kernel $K$, we have the following
joint PDF of the $k$ largest eigenvalues $\lambda_{1}\geq\dots\geq\lambda_{k}$

$f^{(\lambda_{1},\dots,\lambda_{k})}(x_{1},\dots,x_{k})=\det\left(\left[K(x_{i},x_{j})\right]_{i,j=1,\dots,k}\right)\cdot\det\left(I-K^{(x_{1},\dots,x_{k})}\\!\\!\restriction_{(x_{k},\infty)}\right),$

for $x_{1}>\dots>x_{k}$ and vanishes otherwise, where the kernel
$K^{(x_{1},\dots,x_{k})}(x,y)$ is defined as,

$K^{(x_{1},\dots,x_{k})}(x,y)\\!=\\!K(x,y)-\begin{bmatrix}K(x,x_{1})\\\
\vdots\\\
K(x,x_{n})\end{bmatrix}^{T}\\!\\!\begin{bmatrix}K(x_{1},x_{1})&\\!\\!\\!\\!\cdots&\\!\\!\\!\\!K(x_{1},x_{n})\\\
\vdots&\\!\\!\\!\\!\ddots&\\!\\!\\!\\!\vdots\\\
K(x_{n},x_{1})&\\!\\!\\!\\!\cdots&\\!\\!\\!\\!K(x_{n},x_{n})\end{bmatrix}^{-1}\\!\\!\begin{bmatrix}K(x_{1},y)\\\
\vdots\\\ K(x_{n},y)\end{bmatrix}.$

In the following sections we give some examples of applications of Proposition
3.3 with numerical experiments. Furthermore the first row of Figure 1 contains
some visualizations of these eigenvalue statistics that could be obtained from
the joint PDF formula of $k=2$ and $k=3$ extreme eigenvalues.

####  3.4.1. Correlation of the two extreme eigenvalues

In [1] the correlation coefficient $\rho(\lambda_{1},\lambda_{2})$ of the two
largest eigenvalues at the soft-edge is computed, from matrix valued kernels
and generating function approach similar to (3.9). On the other hand, we could
compute $\rho(\lambda_{1},\lambda_{2})$ with the following steps:

  1. (1)

Compute $\mathbb{E}\lambda_{1}\lambda_{2}$ using the joint PDF
$f^{(\lambda_{1},\lambda_{2})}(x_{1},x_{2})$ from Proposition 3.3 and
2-dimensional Gauss quadrature on a (truncated) triangular region.

  2. (2)

Compute $\mathbb{E}\lambda_{1}$, $\mathbb{E}\lambda_{2}$, $\sigma\lambda_{1}$,
$\sigma\lambda_{2}$ using the PDF expressions (1.6), (3.10) and the Gauss-
Legendre quadrature. Infinite intervals such as $(s,\infty)$ are handled by
the strategy explained in Section 7 of [2].

  3. (3)

Compute
$\rho(\lambda_{1},\lambda_{2})=(\mathbb{E}\lambda_{1}\lambda_{2}-\mathbb{E}\lambda_{1}\mathbb{E}\lambda_{2})/(\sigma\lambda_{1}\sigma\lambda_{2})$.

The total computing time for 11 accurate digits,
$\rho(\lambda_{1},\lambda_{2})=0.50564723159$, is 118 seconds. Previously
reported computing time is 16 hours [1]. In addition to the soft-edge, we
compute correlation coefficients of the two smallest eigenvalues at the hard-
edge for $\alpha=0,1,2$. See Table 5 for the result.

| $\rho(\lambda_{1},\lambda_{2})$
---|---
Soft-edge | 0.505 647 231 59
$\alpha=0$ | 0.337 619 085 22
$\alpha=1$ | 0.391 735 693 02
$\alpha=2$ | 0.417 187 915 41
Table 5. Computed values of correlation coefficients of the two largest
eigenvalues at the soft-edge and the two smallest eigenvalues
($\lambda_{1}<\lambda_{2}$) at the hard-edge scaling limit. Computation time
is 117 seconds for the soft-edge and 139 seconds for the whole hard-edge
correlation coefficients. See codes cor-coeff-softedge and cor-coeff-hardedge
for the detailed implementation.

####  3.4.2. First eigenvalue spacing

We compute moments of the distance between the first two eigenvalues (the
first eigenvalue spacing) by computing the PDF and CDF using Proposition 1.1.
In the second row of Figure 1 we also plot some distributions of the first
eigenvalue spacing, using the expressions we derive in this section.

The probability that the first spacing is at least $d$ could be obtained by
integrating over $x$ the probability density of no further eigenvalues
existing in $(x-d,\infty)$ given a level at $x$. From Proposition 1.1 such a
probability density is given as (for example, at the soft-edge)

$K_{\operatorname{Ai}}(x,x)\cdot\det\left(I-K_{\operatorname{Ai}}^{(d)}\\!\\!\restriction_{(x-d,\infty)}\right),$

thus we obtain the CDF $G(s)$ of the first spacing,

$G(d)=1-\int_{\mathbb{R}}K_{\operatorname{Ai}}(x,x)\det(I-K_{\operatorname{Ai}}^{(d)}\\!\\!\restriction_{(x-d,\infty)})dx.$

Moreover, simply using the joint PDF $f^{(\lambda_{1},\lambda_{2})}$ of the
two eigenvalues computed above we obtain

(3.13) $A(d)=\int_{\mathbb{R}}d\cdot f^{(\lambda_{1},\lambda_{2})}(x,x-d)dx,$

which is the PDF of the first spacing. For the implementation we use the
Gauss-Legendre quadrature for the integration.

Mean | Variance | Skewness | Excess Kurtosis
---|---|---|---
1.904 350 489 721 | 0.683 252 055 105 | 0.562 291 976 040 | 0.270 091 960 715
Table 6. The first four moments of the distance between the first two
eigenvalues of the soft-edge scaling, up to 12 digits. Equation (3.13) is used
to compute these values.

Table 6 is the computed first four moments of the first eigenvalue spacing, up
to 12 digits, with a total runtime of 236 seconds. These moments were already
computed in [32] up to $5\sim 9$ digits, with a reported computing time of 5
hours. With (3.13) computation of the moments up to 9 digits needs 29 seconds
of a runtime. Values of $A,G$ of the soft-edge scaling are verified up to 8
digits against Table 2 of [32], with a total computing time of 354 seconds.
Computation of $A,G$ values and comparison with previously known values could
be found in first-spacing.

###  3.5. Sampling Dyson Brownian motion and the Airy process using DPP

In this section we add an additional (time) parameter $t$ as a random matrix
changes through time according to some stochastic processes. A random matrix
diffusion or Dyson process, for example Dyson Browian motion (GUE diffusion),
is another example of DPPs in random matrix theory. A multitime correlation
function for the $N\times N$ GUE diffusion is given in terms of a block matrix
determinant with the extended Hermite kernel $K$ [31],

(3.14) $p(x_{t_{s},i_{s}}\,\,;\,\,s=1,\dots,n\text{ and
}i_{s}=1,\dots,m_{s})=\det\left(\left[K_{j,k}\right]_{j,k=1,\dots,n}\right),$

where $K_{j,k}$ is an $m_{j}\times m_{k}$ matrix,

$K_{j,k}=\left[K(x_{t_{j},i_{j}},x_{t_{k},i_{k}})\right]_{\begin{subarray}{c}i_{j}=1,\dots,m_{j}\\\
i_{k}=1,\dots,m_{k}\end{subarray}}.$

This is essentially the density of eigenvalues of the random matrix stochastic
process at times $\\{t_{s}\\}_{s=1,\dots,n}$ each existing on positions
$x_{t_{s},1},\dots,x_{t_{s},m_{s}}$. The determinantal multitime correlation
function (3.14) also holds for the soft-edge scaling, LUE, and other
orthogonal polynomial ensembles with appropriately computed kernels as proved
in [11]. Indeed the block matrix determinant could be discretized to a block
matrix $K$ kernel for a DPP as prescribed in Section 2.4.

In particular, with the extended Airy kernel

$K_{s,t}^{\text{ext}}(x,y)=\left\\{\begin{array}[]{cc}\int_{0}^{\infty}e^{-\lambda(s-t)}\operatorname{Ai}(x+\lambda)\operatorname{Ai}(y+\lambda)d\lambda&\text{if
}s\geq t\\\
-\int_{-\infty}^{0}e^{-\lambda(s-t)}\operatorname{Ai}(x+\lambda)\operatorname{Ai}(y+\lambda)d\lambda&\text{if
}s<t\end{array}\right.,$

one has the multitime correlation function of the soft-edge scaling limit by
(3.14). A special interest lies in the largest eigenvalue of this process and
is called the Airy2 process, or just simply, the Airy process. Obviously
extending (2.2), the largest eigenvalue process is described by the following
Fredholm determinant,

(3.15) $\mathbb{P}(\text{$\mathcal{A}(t_{1})\leq
s_{1},\dots,\mathcal{A}(t_{n})\leq
s_{n}$})=\det\left(I-K\\!\\!\restriction_{L^{2}(s_{1},\infty)\oplus\cdots\oplus
L^{2}(s_{n},\infty)}\right),$

where $K$ is the block kernel given as

$K=\begin{bmatrix}K_{t_{1},t_{1}}^{\text{ext}}&\cdots&K_{t_{1},t_{n}}^{\text{ext}}\\\
\vdots&\ddots&\vdots\\\
K_{t_{n},t_{1}}^{\text{ext}}&\cdots&K_{t_{n},t_{n}}^{\text{ext}}\end{bmatrix}.$

The Airy process is related to a number of applications, including the
polynuclear growth process [19, 28], the NPR boundary process of the Aztec
diamond domino tiling [20] (i.e., the top DR path as $n\to\infty$ in Section
2.3), totally asymmetric simple exclusion process (TASEP), corner growth
process [19], and eventually related to the KPZ universality class with the
narrow wedge initial condition.

An example of numerical experiments on the Airy process is [2], where
Bornemann use the $2\times 2$ block Fredholm determinant (3.15) to compute the
two-point covariance $\text{cov}(\mathcal{A}(0),\mathcal{A}(t))$ and compare
those values against its large $t$ and small $t$ expansions.

Figure 5. Sketch of the sampling technique used to sample the Airy process.
Each block represents discretized $K_{t_{i},t_{j}}^{\text{ext}}$. At each
timestep $t_{i}$ we perform Algorithm 3, starting from $+\infty$, until we hit
the first (largest) eigenvalue. Then we jump to the next timestep,
disregarding (not observing) rest of the eigenvalues at the current timestep.
Note that each diagonal block of the kernel is the usual Airy kernel (3.2).

Here we demonstrate another numerical experiment on the Airy process: Sampling
Airy processes through multitime correlation function (3.14) and the DPP
sampler. We use $201\times 201$ block matrix for sampling. The extended Airy
kernel is non-Hermitian and non-Projection, which could only be sampled by
Algorithm 3. A simple multi-step modification of Algorithm 3 is handy; In each
timestep $t_{i}$, we proceed from $+\infty$ to $-\infty$, and jump to next
block (next timestep $t_{i+1}$) when we find a largest eigenvalue at the
current timestep. See Figure 5 for the illustration. More details such as
discretization and truncating the interval to get finite kernel could be found
in Section 2.4.

In the left of Figure 6 is the plot of five samples of the Airy process,
sampled from the DPP defined by multitime correlation function of the (soft-
edge scaling) GUE diffusion, i.e., the extended Airy kernel. These are exact
in a sense that the samples are not large $N$ asypmtotics where $N\to\infty$
is the Airy process. Such large $N$ approximations are drawn in red in the
right of Figure 6.

Figure 6. Five samples of the Airy process from $t=0$ to $t=5$, sampled at
interval points with $dt=0.025$. We truncated eigenvalue space to
$[-5.0,2.5]$, as probabilities that a sample from the Tracy–Widom distribution
larger than 2.5 and smaller than -5 are both around $2\times 10^{-5}$.
Eigenvalue domain is then discretized into 150 points, finally yielding
$30150\times 30150$ kernel for the DPP. Sampling time for a single sample is
around 4 hours. On the right we additionally draw five samples of the largest
eigenvalue process $\lambda_{\text{max}}(t)$ of the Dyson Brownian motion with
$N=200$, recentered and rescaled according to
$\sqrt{2}N^{\frac{1}{6}}\left(\lambda_{\text{max}}(N^{-\frac{1}{3}}t)-\sqrt{2N}\right),\hskip
56.9055pt$ where it is known that as $N\to\infty$ the above scaling converges
to the Airy process.

### Numerical experiment details

All codes mentioned in the paper could be found online:
https://github.com/sw2030/RMTexperiments/codes. For numerical experiments
(except Section 3.5) discussed in this work, we used a single core of an Apple
M1 Pro CPU. For the Airy process sampling discussed in Section 3.5, we used 64
cores from four Xeon P8 CPUs for computing the DPP kernel and a single core of
Xeon P8 CPU for sampling, from the MIT Supercloud server [29].

### Acknowledgements

This material is based upon work supported by the National Science Foundation
under Grant No. DMS-1926686. The authors acknowledge the MIT SuperCloud and
Lincoln Laboratory Supercomputing Center for providing HPC resources that have
contributed to the research results reported within this paper. This material
is based upon work supported by the National Science Foundation under grant
no. OAC-1835443, grant no. SII-2029670, grant no. ECCS-2029670, grant no.
OAC-2103804, and grant no. PHY-2021825. We also gratefully acknowledge the
U.S. Agency for International Development through Penn State for grant no.
S002283-USAID. The information, data, or work presented herein was funded in
part by the Advanced Research Projects Agency-Energy (ARPA-E), U.S. Department
of Energy, under Award Number DE-AR0001211 and DE-AR0001222. We also
gratefully acknowledge the U.S. Agency for International Development through
Penn State for grant no. S002283-USAID. The views and opinions of authors
expressed herein do not necessarily state or reflect those of the United
States Government or any agency thereof. This material was supported by The
Research Council of Norway and Equinor ASA through Research Council project
“308817 - Digital wells for optimal production and drainage”. Research was
sponsored by the United States Air Force Research Laboratory and the United
States Air Force Artificial Intelligence Accelerator and was accomplished
under Cooperative Agreement Number FA8750-19-2-1000. The views and conclusions
contained in this document are those of the authors and should not be
interpreted as representing the official policies, either expressed or
implied, of the United States Air Force or the U.S. Government. The U.S.
Government is authorized to reproduce and distribute reprints for Government
purposes notwithstanding any copyright notation herein.

## References

  * [1] Folkmar Bornemann.  On the numerical evaluation of distributions in random matrix theory: a review.  Markov Process. Relat. Fields, 16(4):803–866, 2010.
  * [2] Folkmar Bornemann.  On the numerical evaluation of Fredholm determinants.  Mathematics of Computation, 79(270):871–915, 2010.
  * [3] Folkmar Bornemann.  A Stirling-type formula for the distribution of the length of longest increasing subsequences.  Foundations of Computational Mathematics, pages 1–39, 2023.
  * [4] Alexei Borodin and Peter J Forrester.  Increasing subsequences and the hard-to-soft edge transition in matrix ensembles.  Journal of Physics A: Mathematical and General, 36(12):2963, 2003\.
  * [5] Alexei Borodin and Eric M Rains.  Eynard–Mehta theorem, Schur process, and their Pfaffian analogs.  Journal of statistical physics, 121:291–317, 2005.
  * [6] Sunil Chhita, Kurt Johansson, and Benjamin Young.  Asymptotic domino statistics in the Aztec diamond.  The Annals of Applied Probability, 25(3):1232–1278, 2015.
  * [7] Michal Derezinski, Daniele Calandriello, and Michal Valko.  Exact sampling of determinantal point processes with sublinear time preprocessing.  Advances in neural information processing systems, 32, 2019.
  * [8] Freeman J Dyson.  Statistical theory of the energy levels of complex systems III.  Journal of Mathematical Physics, 3(1):166–175, 1962.
  * [9] Alan Edelman.  MIT class 18.338 Eigenvalues of random matrices, class notes, Fall 2022.
  * [10] Noam Elkies, Greg Kuperberg, Michael Larsen, and James Propp.  Alternating-sign matrices and domino tilings (Part II).  Journal of Algebraic Combinatorics, 1:219–234, 1992.
  * [11] Bertrand Eynard and Madan Lal Mehta.  Matrices coupled in a chain: I. eigenvalue correlations.  Journal of Physics A: Mathematical and General, 31(19):4449, 1998\.
  * [12] Peter J Forrester.  The spectrum edge of random matrix ensembles.  Nuclear Physics B, 402(3):709–728, 1993.
  * [13] Peter J Forrester.  Log-Gases and Random Matrices (LMS-34).  Princeton University Press, 2010.
  * [14] Peter J Forrester and Nicholas S Witte.  The distribution of the first eigenvalue spacing at the hard edge of the Laguerre unitary ensemble.  Kyushu Journal of Mathematics, 61(2):457–526, 2007.
  * [15] Michel Guillaume and Jean Gautier.  On sampling determinantal point processes.  PhD thesis, Centrale Lille Institut, 2020.
  * [16] J Ben Hough, Manjunath Krishnapur, Yuval Peres, and Bálint Virág.  Determinantal processes and independence.  Probability surveys, 3:206–229, 2006.
  * [17] William Jockusch, James Propp, and Peter Shor.  Random domino tilings and the arctic circle theorem.  arXiv preprint math/9801068, 1998.
  * [18] Kurt Johansson.  Non-intersecting paths, random tilings and random matrices.  Probability theory and related fields, 123(2):225–280, 2002.
  * [19] Kurt Johansson.  Discrete polynuclear growth and determinantal processes.  Communications in Mathematical Physics, 242:277–329, 2003.
  * [20] Kurt Johansson.  The arctic circle boundary and the Airy process.  The Annals of Probability, 33(1):1–30, 2005.
  * [21] Richard Kenyon.  Local statistics of lattice dimers.  In Annales de l’Institut Henri Poincare (B) Probability and Statistics, volume 33, pages 591–618. Elsevier, 1997.
  * [22] Alex Kulesza and Ben Taskar.  k-DPPs: Fixed-size determinantal point processes.  In Proceedings of the 28th International Conference on Machine Learning (ICML-11), pages 1193–1200, 2011.
  * [23] Alex Kulesza and Ben Taskar.  Learning determinantal point processes.  In Proceedings of the Twenty-Seventh Conference on Uncertainty in Artificial Intelligence, pages 419–427, 2011.
  * [24] Claire Launay, Bruno Galerne, and Agnès Desolneux.  Exact sampling of determinantal point processes without eigendecomposition.  Journal of Applied Probability, 57(4):1198–1221, 2020.
  * [25] Madan Lal Mehta.  Random Matrices.  Elsevier, 2004.
  * [26] Taro Nagao and Miki Wadati.  Correlation functions of random matrix ensembles related to classical orthogonal polynomials.  Journal of the Physical Society of Japan, 60(10):3298–3322, 1991\.
  * [27] Jack Poulson.  High-performance sampling of generic determinantal point processes.  Philosophical Transactions of the Royal Society A, 378(2166):20190059, 2020.
  * [28] Michael Prähofer and Herbert Spohn.  Scale invariance of the PNG droplet and the Airy process.  Journal of statistical physics, 108:1071–1106, 2002.
  * [29] Albert Reuther, Jeremy Kepner, Chansup Byun, Siddharth Samsi, William Arcand, David Bestor, Bill Bergeron, Vijay Gadepally, Michael Houle, Matthew Hubbell, Michael Jones, Anna Klein, Lauren Milechin, Julia Mullen, Andrew Prout, Antonio Rosa, Charles Yee, and Peter Michaleas.  Interactive supercomputing on 40,000 cores for machine learning and data analysis.  In 2018 IEEE High Performance extreme Computing Conference (HPEC), pages 1–6. IEEE, 2018.
  * [30] Richard P. Stanley.  Enumerative Combinatorics, volume 2.  Cambridge University Press, 1999.
  * [31] Craig A Tracy and Harold Widom.  Differential equations for Dyson processes.  Communications in mathematical physics, 252:7–41, 2004.
  * [32] Nicholas S Witte, Folkmar Bornemann, and Peter J Forrester.  Joint distribution of the first and second eigenvalues at the soft edge of unitary ensembles.  Nonlinearity, 26(6):1799, 2013.