RandomMatrices.jl
Random matrix package for Julia.
This extends the Distributions package to provide methods for working with matrixvalued random variables, a.k.a. random matrices. State of the art methods for computing random matrix samples and their associated distributions are provided.
The names of the various ensembles can vary widely across disciplines. Where possible, synonyms are listed.
Additional functionality is provided when these optional packages are installed:
 Symbolic manipulation of Haar matrices with GSL.jl
Gaussian matrix ensembles
Much of classical random matrix theory has focused on matrices with matrix elements comprised of
independently and identically distributed (iid) real, complex or quaternionic Gaussians.
(Traditionally, these are associated with a parameter beta
tracking the number of independent
real random variables per matrix element, i.e. beta=1,2,4
respectively. This is also referred
to as the Dyson 3fold way.)
Methods are provided for calculating random variates (samples) and various properties of these
random matrices.
The hierarchy of dense matrices provided are
 Ginibre ensemble  all matrix elements are iid with no global symmetry
 Hermite ensemble  one global symmetry
 Gaussian orthogonal ensemble (GOE,
beta=1
)  real and symmetric  Gaussian unitary ensemble (GUE,
beta=2
)  complex and Hermitian  Gaussian symplectic ensemble (GSE,
beta=4
)  quaternionic and selfdual
 Gaussian orthogonal ensemble (GOE,
 Circular ensemble  uniformly distributed with
det=1
 Circular orthogonal ensemble (COE,
beta=1
)  Circular unitary ensemble (CUE,
beta=2
)  Circular symplectic ensemble (CSE,
beta=4
)
 Circular orthogonal ensemble (COE,
 Laguerre matrices = white Wishart matrices
 Jacobi matrices = MANOVA matrices
Unless otherwise specified, beta=1,2,4
are supported. For the symplectic matrices beta=4
,
the 2x2 outer blockdiagonal complex representation USp(2N)
is used.
Joint probability density functions (jpdfs)
Given eigenvalues lambda
and the beta
parameter of the random matrix distribution:
VandermondeDeterminant(lambda, beta)
computes the Vandermonde determinantHermiteJPDF(lambda, beta)
computes the jpdf for the Hermite ensembleLaguerreJPDF(lambda, n, beta)
computes the jpdf for the Laguerre(n) ensembleJacobiJPDF(lambda, n1, n2, beta)
computes the jpdf for the Jacobi(n1, n2) ensemble
Matrix samples
Constructs samples of random matrices corresponding to the classical Gaussian Hermite, Laguerre(m) and Jacobi(m1, m2) ensembles.

GaussianHermiteMatrix(n, beta)
,GaussianLaguerreMatrix(n, m, beta)
,GaussianJacobiMatrix(n, m1, m2, beta)
each construct a sample densen
xn
matrix for the corresponding matrix ensemble withbeta=1,2,4

GaussianHermiteTridiagonalMatrix(n, beta)
,GaussianLaguerreTridiagonalMatrix(n, m, beta)
,GaussianJacobiSparseMatrix(n, m1, m2, beta)
each construct a sparsen
xn
matrix for the corresponding matrix ensemble for arbitrary positive finitebeta
.GaussianHermiteTridiagonalMatrix(n, Inf)
is also allowed. These sampled matrices have the same eigenvalues as above but are much faster to diagonalize oweing to their sparsity. They also extend Dyson's threefold way to arbitrarybeta
. 
GaussianHermiteSamples(n, beta)
,GaussianLaguerreSamples(n, m, beta)
,GaussianJacobiSamples(n, m1, m2, beta)
return a set ofn
eigenvalues from the sparse random matrix samples 
HaarMatrix(n, beta)
Generates a random matrix from thebeta
circular ensemble. 
HaarMatrix(n, beta, correction)
provides finegrained control of what kind of correction is applied to the raw QR decomposition. By default,correction=1
(Edelman's correction) is used. Other valid values are0
(no correction) and2
(Mezzadri's correction). 
NeedsPiecewiseCorrection()
implements a simple test to see if a correction is necessary.
The parameters m
, m1
, m2
refer to the number to independent "data" degrees of freedom.
For the dense samples these must be Integer
s but can be Real
s for the rest.
Formal power series
Allows for manipulations of formal power series (fps) and formal Laurent series (fLs), which come in handy for the computation of free cumulants.
Types
FormalPowerSeries
: power series with coefficients allowed only for nonnegative integer powersFormalLaurentSeries
: power series with coefficients allowed for all integer powers
FormalPowerSeries methods
Elementary operations
 basic arithmetic operations
==
,+
,
,^
*
computes the Cauchy product (discrete convolution).*
computes the Hadamard product (elementwise multiplication)compose(P,Q)
computes the series composition P.Qderivative
computes the series derivativereciprocal
computes the series reciprocal
Utility methods
trim(P)
removes extraneous zeroes in the internal representation ofP
isalmostunit(P)
determines ifP
is an almost unit seriesisconstant(P)
determines ifP
is a constant seriesisnonunit(P)
determines ifP
is a nonunit seriesisunit(P)
determines ifP
is a unit seriesMatrixForm(P)
returns a matrix representation ofP
as an upper triangular Toeplitz matrixtovector
returns the series coefficients
Densities
Famous distributions in random matrix theory
Semicircle
provides the semicircle distributionTracyWidom
computes the TracyWidom density distribution by bruteforce integration of the Painlevé II equation
Utility functions
hist_eig
computes the histogram of eigenvalues of a matrix using the method of Sturm sequences. This is recommended forSymTridiagonal
matrices as it is significantly faster thanhist(eigvals())
This is also implemented for dense matrices, but it is pretty slow and not really practical.
Stochastic processes
Julia iterators for stochastic operators.
All subtypes of StochasticProcess
contain at least one field, dt
,
representing the time interval being discretized over.
The available StochasticProcess
es are
BrownianProcess(dt)
: Brownian random walk. The state of the iterator is the cumulative displacement of the random walk.WhiteNoiseProcess(dt)
: White noise. The value of this iterator israndn()*dt
. The state associated with this iterator isnothing
.StochasticAiryProcess(dt, beta)
: stochastic Airy process with real positivebeta
. The value of this iterator in the limit of an infinite number of iterations is known to follow thebeta
TracyWidom law. The state associated with this iteratior is aSymTridiagonal
matrix whose largest eigenvalue is the value of this process.
Invariant ensembles
InvariantEnsemble(str,n)
supports n x n unitary invariant ensemble
with distribution. This has been moved to separate package InvariantEnsembles.jl
References

James Albrecht, Cy Chan, and Alan Edelman, "Sturm Sequences and Random Eigenvalue Distributions", Foundations of Computational Mathematics, vol. 9 iss. 4 (2009), pp 461483. [pdf] [doi]

Ioana Dumitriu and Alan Edelman, "Matrix Models for Beta Ensembles", Journal of Mathematical Physics, vol. 43 no. 11 (2002), pp. 58305547 [[doi]](http://dx.doi.org/doi: 10.1063/1.1507823) arXiv:mathph/0206043

Alan Edelman, PerOlof Persson and Brian D Sutton, "The fourfold way", Journal of Mathematical Physics, submitted (2013). [pdf]

Alan Edelman and Brian D. Sutton, "The betaJacobi matrix model, the CS decomposition, and generalized singular value problems", Foundations of Computational Mathematics, vol. 8 iss. 2 (2008), pp 259285. [pdf] [doi]

Peter Henrici, Applied and Computational Complex Analysis, Volume I: Power SeriesIntegrationConformal MappingLocation of Zeros, WileyInterscience: New York, 1974 [worldcat]

Frank Mezzadri, "How to generate random matrices from the classical compact groups", Notices of the AMS, vol. 54 (2007), pp592604 [arXiv]