8 days ago
Random matrix package for Julia.
This extends the Distributions package to provide methods for working with matrix-valued 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:
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 3-fold way.)
Methods are provided for calculating random variates (samples) and various properties of these
The hierarchy of dense matrices provided are
beta=1) - real and symmetric
beta=2) - complex and Hermitian
beta=4) - quaternionic and self-dual
Unless otherwise specified,
beta=1,2,4 are supported. For the symplectic matrices
the 2x2 outer block-diagonal complex representation
USp(2N) is used.
lambda and the
beta parameter of the random matrix distribution:
VandermondeDeterminant(lambda, beta)computes the Vandermonde determinant
HermiteJPDF(lambda, beta)computes the jpdf for the Hermite ensemble
LaguerreJPDF(lambda, n, beta)computes the jpdf for the Laguerre(n) ensemble
JacobiJPDF(lambda, n1, n2, beta)computes the jpdf for the Jacobi(n1, n2) ensemble
Constructs samples of random matrices corresponding to the classical Gaussian Hermite, Laguerre(m) and Jacobi(m1, m2) ensembles.
GaussianLaguerreMatrix(n, m, beta),
GaussianJacobiMatrix(n, m1, m2, beta)
each construct a sample dense
n matrix for the corresponding matrix
GaussianLaguerreTridiagonalMatrix(n, m, beta),
GaussianJacobiSparseMatrix(n, m1, m2, beta)
each construct a sparse
n matrix for the corresponding matrix ensemble
for arbitrary positive finite
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 arbitrary
GaussianLaguerreSamples(n, m, beta),
GaussianJacobiSamples(n, m1, m2, beta)
return a set of
n eigenvalues from the sparse random matrix samples
Generates a random matrix from the
HaarMatrix(n, beta, correction)provides fine-grained 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 are
0(no correction) and
NeedsPiecewiseCorrection()implements a simple test to see if a correction is necessary.
Generates a unitary invariant ensemble, where str determines the
potential of the ensemble, see below.
Only available if ApproxFun package is installed.
m2 refer to the number to independent "data" degrees of freedom.
For the dense samples these must be
Integers but can be
Reals for the rest.
Allows for manipulations of formal power series (fps) and formal Laurent series (fLs), which come in handy for the computation of free cumulants.
FormalPowerSeries: power series with coefficients allowed only for non-negative integer powers
FormalLaurentSeries: power series with coefficients allowed for all integer powers
*computes the Cauchy product (discrete convolution)
.*computes the Hadamard product (elementwise multiplication)
compose(P,Q)computes the series composition P.Q
derivativecomputes the series derivative
reciprocalcomputes the series reciprocal
trim(P)removes extraneous zeroes in the internal representation of
Pis an almost unit series
Pis a constant series
Pis a non-unit series
Pis a unit series
MatrixForm(P)returns a matrix representation of
Pas an upper triangular Toeplitz matrix
tovectorreturns the series coefficients
Famous distributions in random matrix theory
Semicircleprovides the semicircle distribution
TracyWidomcomputes the Tracy-Widom density distribution by brute-force integration of the Painlevé II equation
hist_eigcomputes the histogram of eigenvalues of a matrix using the method of Sturm sequences. This is recommended for
SymTridiagonalmatrices as it is significantly faster than
hist(eigvals())This is also implemented for dense matrices, but it is pretty slow and not really practical.
Julia iterators for stochastic operators.
All subtypes of
StochasticProcess contain at least one field,
representing the time interval being discretized over.
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 is
randn()*dt. The state associated with this iterator is
StochasticAiryProcess(dt, beta): stochastic Airy process with real positive
beta. The value of this iterator in the limit of an infinite number of iterations is known to follow the
beta-Tracy-Widom law. The state associated with this iteratior is a
SymTridiagonalmatrix whose largest eigenvalue is the value of this process.
InvariantEnsemble(str,n) supports n x n unitary invariant ensemble
exp(- Tr Q(M)) dM
str specifies an ensemble with precomputed recurrence coefficients.
The currently include ensembles are
||n (M^4/20 - 4/15M^3 +M^2/5 + 8/5M)|
Alan Edelman, Per-Olof Persson and Brian D Sutton, "The fourfold way", Journal of Mathematical Physics, submitted (2013). [pdf]
Alan Edelman and Brian D. Sutton, "The beta-Jacobi matrix model, the CS decomposition, and generalized singular value problems", Foundations of Computational Mathematics, vol. 8 iss. 2 (2008), pp 259-285. [pdf] [doi]
Peter Henrici, Applied and Computational Complex Analysis, Volume I: Power Series---Integration---Conformal Mapping---Location of Zeros, Wiley-Interscience: New York, 1974 [worldcat]
Frank Mezzadri, "How to generate random matrices from the classical compact groups", Notices of the AMS, vol. 54 (2007), pp592-604 [arXiv]