1
18.338J/16.394J: The Mathematics of In?nite Random Matrices
Experiments with the Classical Ensembles
Professor Alan Edelman
Handout #3, Tuesday, September 14, 2004
The Classical Random Matrix Ensembles
The Wigner Matrix (or Hermite Ensemble)
The Wigner matrices [1, 2] are often known as the Hermite or Gaussian ensembles are well studied in
physics and in the book by Mehta [3]. The term Wigner matrix does not require the entries be normal,
though the term Gaussian ensemble and Hermite ensemble does. Typically what is required is that the
elements have the same independence properties and variance as the Gaussian case. This means variance
1/n o? the diagonal and sqrt(2)/n on the diagonal, and elements independent on the upper triangle.
Let G be a N ×N random matrix with independent, zero mean, unit variance elements. The Wigner
matrix W can be obtained from G as simply
′
G + G
W = . (1) √
2N
Code 1 lists a MATLAB function for producting a Wigner matrix from i.i.d. standard normals (randn in
MATLAB notation).
The speci?c names Gaussian orthogonal ensemble, Gaussian unitary ensemble, and Gaussian symplectic
ensemble or GOE, GUE, GSE were named by Dyson [4] and refer to the real, complex, and quaternion
cases respectively. These names refer to invariance properties of the distribution.
The name Hermite ensemble refers to the connection to the Hermite weight function and the associated
Hermite orthogonal polynomials that are so intimitely connected to the ?nite theory.
? ?
? ?
function w = wigner(n,isreal);
%WIGNER The Wigner matrix.
% WIGNER(N,ISREAL) generates an N x N symettric Wigner matrix.
% If ISREAL = 1 then the elements of W are real.
% If ISREAL = 0 then the elements of W are complex.
%
% N is an integer while ISREAL equals either 0 or 1.
% W is an N x N matrix.
%
% References:
% [1] Alan Edelman, Handout 3: Experiments with Classical
% Matrix Ensembles, Fall 2004,
% Course Notes 18.338.
% [2] Alan Edelman, Random Matrix Eigenvalues.
% [3] E. P. Wigner, Characteristic vectors of bordered matrices with
% infinite dimensions, Annals of Mathematics,
% vol. 62, pages 548-564, 1955.
%
% Alan Edelman and Raj Rao, Sept. 2004.
% $Revision: 1.0 $ $Date: 2004/09/10 23:21:18 $
if(isreal==1)
g = randn(n,n);
else
g = (randn(n,n) + i*randn(n,n))/sqrt(2);
end
w = (g + g’) / sqrt(2*n);
Code 1
The Wishart Matrix (or Laguerre Ensemble)
The Wishart matrices [5], also sometimes referred to as Grammian matrices [6] or sample covariance
matrices, are well studied in statistics and in the book by Muirhead [7].
Let G be a N ×M random matrix with independent, zero mean, unit variance elements. The Wishart
matrix W can be obtained from G as simply
1 ′
W = GG . (2)
M
The term Wishart matrix does require that the entries of G are i.i.d. normally distributed elements as does
the term Laguerre ensemble. We shall refer to this as a “pure” Wishart matrix when we need to make that
distinction. When analyzing the limiting density, typically what is required is that the elements of G are
i.i.d. with zero mean and unit variance. Associated with the Wishart matrix is the parameter c = N/M
which is simply the ratio of the number of rows to the number of columns of the matrix G. We shall
employ the notation W(c) to denote the generalized Wishart matrix with its associated parameter c.
Code 2 lists a MATLAB function for producing a “pure” Wishart matrix from i.i.d. standard normals (randn
in MATLAB notation).
The name Laguerre ensemble refers to the connection to the Laguerre weight function and the associated
Laguerre orthogonal polynomials that are so intimately connected to the ?nite theory.
? ?
? ?
function w = wishart(n,m,isreal);
%WISHART The Wishart matrix.
% WISHART(N,M,ISREAL) generates an N x N Hermitian Wishart matrix.
% If ISREAL = 1 then the elements of W are real.
% If ISREAL = 0 then the elements of W are complex.
%
% N and M are integers while ISREAL equals either 0 or 1.
% W is an N x N matrix.
%
% References:
% [1] Alan Edelman, Handout 3: Experiments with Classical
% Matrix Ensembles, Fall 2004,
% Course Notes 18.338.
% [2] Alan Edelman, Random Matrix Eigenvalues.
% [3] J. Wishart, The generalized product moment distribution in
% samples from a normal multivariate population,
% Biometrika, vol. 20-A, pages 32-52, 1928.
%
% Alan Edelman and Raj Rao, Sept. 2004.
% $Revision: 1.0 $ $Date: 2004/09/10 23:45:18 $
if(isreal==1)
g = randn(n,m);
else
g = (randn(n,m) + i*randn(n,m))/sqrt(2);
end
w = (g * g’) / m;
Code 2
The MANOVA Matrix (or Jacobi Ensemble)
The name Jacobi ensemble refers to the connection to the Jacobi weight function and the associated Jacobi
orthogonal polynomials that are so intimately connected to the ?nite theory.
The MANOVA matrix, J, can be de?ned in terms of two Wishart matrices, W(c
1
) and W(c
2
), as simply
J = (W(c
1
) + W(c
2
))
?1
W(c
1
). (3)
The term Jacobi ensemble requires that the Wishart matrices, W(c
1
) and W(c
2
), be “pure”. We shall
employ the notation J(c
1
, c
2
) to denote the MANOVA matrix where c
1
= N/M
1
and c
2
= N/M
2
are the
associated parameters of the two Wishart matrices, respectively. Note that the de?nition of the MANOVA
matrix implies that c
1
= N/M
1
< 1 i.e. M
1
> N.
? ?
? ?
function j = manova(n,m1,m2,isreal);
%MANOVA The MANOVA matrix.
% MANOVA(N,M1,M2,ISREAL) generates an N x N MANOVA matrix.
% If ISREAL = 1 then the elements of W are real.
% If ISREAL = 0 then the elements of W are complex.
%
% N , M1 , and M2 are integers while ISREAL equals either 0 or 1.
% N < M1.
% J is an N x N matrix.
%
% References:
% [1] Alan Edelman, Handout 3: Experiments with Classical
% Matrix Ensembles, Fall 2004,
% Course Notes 18.338.
% [2] Alan Edelman, Random Matrix Eigenvalues.
% [3] R. J. Muirhead, Aspects of Multivariate Statistical Theory,
% John Wiley & Sons, New York, 1982.
%
%
% Alan Edelman and Raj Rao, Sept. 2004.
% $Revision: 1.0 $ $Date: 2004/09/10 23:55:18 $
w1 = wishart(n,m1,isreal); % Generate Wishart Matrices
w2 = wishart(n,m2,isreal);
j = (w1 + w2) w1;
Code 3
Code 3 lists a MATLAB function for producing a MANOVA matrix from pure Wishart matrices.
? ?
?
? ?
?
? ?
? ?
? ?
2 Limiting densities
The limiting density of the eigenvalues of a Wigner matrix [1, 2] is given by
√
4 ?x
2
Density = . (4)
2π
The limiting density of the eigenvalues of a Wishart matrix [8] is given by
1
Density = max 0, 1 ? δ(x) +
(x ?b
?
)(b
+
?x)
I
[ b
?
,b
+
]
(5)
c 2πxc
where b
±
= (1 ±
√
c)
2
and I is the indicator function.
The limiting density of the eigenvalues of a MANOVA matrix is given by
Density = max 0, 1 ?
1
δ(x ?1) +
(x ?b
?
)(b
+
?x)
I
[ b
?
,b
+
]
(6)
c
2
2πc
0
(x)
2 3
where c
0
(x) = c
1
x + x
3
c
1
?2c
1
x ?c
2
x + c
2
x
2
and b
±
is given by
b
±
=
2
(1 ?c
1
)
2
. (7)
c
1
?c1 + 2 + c
2
?c
1
c
2
?2
√
c
1
+ c
2
?c
1
c
2
Verify that the MATLAB functions listed do indeed compute these densities symbolically.
function f = wigneredf;
%WIGNEREDF Density of eigenvalues of an infinite Wigner matrix.
% WIGNEREDF returns the theoretical density of the eigenvalues of
% an infinite Wigner matrix
%
% F is a symbolic variable.
%
% References:
% [1] Alan Edelman, Handout 3: Experiments with Classical
% Matrix Ensembles, Fall 2004,
% Course Notes 18.338.
% [2] Alan Edelman, Random Matrix Eigenvalues.
% [3] E. P. Wigner, Characteristic vectors of bordered matrices with
% infinite dimensions, Annals of Mathematics,
% vol. 62, pages 548-564, 1955.
%
% Alan Edelman and Raj Rao, Sept. 2004.
% $Revision: 1.0 $ $Date: 2004/09/11 20:47:15 $
syms f x
f = sqrt(4-x^2)/(2*pi);
Code 4
? ?
? ?
function f = wishartedf(c)
%WISHARTEDF Density of eigenvalues of an infinite Wishart matrix.
% THWISHART(C) returns the theoretical limiting density of the
% eigenvalues of an infinite Wishart matrix with
% C = N / M where N and M are the parameters of
% the Wishart matrix.
%
% C is a positive real number less than 1.
% F is a symbolic variable.
%
% References:
% [1] Alan Edelman, Handout 3: Experiments with Classical
% Matrix Ensembles, Fall 2004,
% Course Notes 18.338.
% [2] Alan Edelman, Random Matrix Eigenvalues.
% [3] J. Wishart, The generalized product moment distribution in
% samples from a normal multivariate population,
% Biometrika, vol. 20-A, pages 32-52, 1928.
%
% Alan Edelman and Raj Rao, Sept. 2004.
% $Revision: 1.0 $ $Date: 2004/09/11 20:49:20 $
syms f x
a1 = (1-sqrt(c))^2;
a2 = (1+sqrt(c))^2;
f = sqrt((x-a1)*(a2-x))/(2*pi*x*c);
Code 5
? ?
? ?
3
function f = manovaedf(c1,c2);
%MANOVAEDF Theoretical density of eigenvalues of an infinite MANOVA matrix.
% MANOVAEDF(C1,C2) returns the limiting theoretical density of the
% eigenvalues of an infinite MANOVA matrix,
% where C1 = N / M1 and C2 = N / M2 and N1, M1, and
% M2 are parameters of the MANOVA matrix.
%
% C1 and C2 are positive real numbers less than 1.
% F is a symbolic variable.
%
% References:
% [1] Alan Edelman, Handout 3: Experiments with Classical
% Matrix Ensembles, Fall 2004,
% Course Notes 18.338.
% [2] Alan Edelman, Random Matrix Eigenvalues.
% [3] R. J. Muirhead, Aspects of Multivariate Statistical Theory,
% John Wiley & Sons, New York, 1982.
%
%
% Alan Edelman and Raj Rao, Sept. 2004.
% $Revision: 1.0 $ $Date: 2004/09/10 23:55:18 $
syms x f
b0 = c1*x-c2*x-c1+2;
b1 = -2*c2*x^2+2*x-3*c1*x+c1+c2*x-1+2*c1*x^2;
b2 = c1*x-2*c1*x^2+c2*x^2-x^3*c2+x^3*c1;
f = sqrt(4*b2*b0-b1^2)/(2*pi*b2);
Code 6
Experiments
.
Put all
in the same directory and run the script classicexpt. Try it with di?erent sizes of
matrices and see how quickly the histograms converge.
code
? ?
? ?
%classicexpt.m
%Code 1.2 of Random Eigenvalues by Alan Edelman and Raj Rao
%Experiment: Generate many realizations of the classical random matrices.
%Observation: Histogram the eigenvalues of the random matrices.
%Theory: Falls on the predicted curves.
%Reference: [1] Alan Edelman, Handout 3, Experiments with Classical
% Matrix Ensembles, Fall 2004.
% [2] Alan Edelman, Random Matrix Eigenvalues.
trials = 1000; % Number of trials
e = []; % Array for collecting eigenvalues
syms f x; % Symbolic variables for theoretical predictions
% Requires Symbolic Toolbox in MATLAB
figure; set(gcf,’DoubleBuffer’,’On’);
warning off;
for i = 1 : trials
% Generate random matrix (uncomment the appropriate line)
isreal = 1;
% N=50; C = wigner(N,isreal); f = wigneredf;
% N=50; M=100; C = wishart(N,M,isreal); f = wishartedf(N/M);
% N=50; M1=100; M2=100; C = manova(N,M1,M2,isreal); f = manovaedf(N/M1,N/M2);
% Collect the eigenvalues
e = [e;real(eig(C))]; % Take real part to compensate for numerical issues
% Plot density of eigenvalues
x0=-2; binsize=0.05; xf=4;
[h,hn,xspan] = histn(e,x0,binsize,xf);
axis([x0 xf 0 1.5]); hold on
% Plot theory
F = real(subs(f,xspan)); F(F==inf)=0;
hold on
ts=[’# trials = ’ num2str(i)]; title(ts,’FontSize’,15);
plot(xspan,F,’red’,’LineWidth’,2);drawnow;
hold off
end
Code 7
4 Exercises
Try the following by appropriately modifying the scripts/functions described in this handout:
? Do the Gaussian matrices have to be circular for the limiting density expressions to be valid? When
isreal = 0 let the real and the complex parts have di?erent variances (keep the sum of the variances
equal to 1!). What happens?
? Try other distributions such as sign(randn), the exponential distribution. What happens?
? Try a distribution with non-zero mean. What can you say about the limiting density?
? Do the number of trials or the size of the matrix at which the density “starts converging” depend on
the distribution? Or on whether it is real or complex?
? Compute the ?rst k moments of these ensembles. Plot the mean and the variance of these moments
as a function of the matrix size N.
References
[1] Eugene P. Wigner. Characteristic vectors of bordered matrices with in?nite dimensions. Annals of
Math., 62:548–564, 1955.
[2] Eugene P. Wigner. On the distribution of the roots of certain symmetric matrices. Annals of
Mathematics, 67:325–327, 1958.
[3] Madan Lal Mehta. Random Matrices. Academic Press, Boston, second edition, 1991.
[4] Freeman J. Dyson. A class of matrix ensembles. J. Math. Phys., 13:90–97, 1972.
[5] J. Wishart. The generalized product moment distribution in samples from a normal multivariate
population. Biometrika, 20 A:32–52, 1928.
[6] V. L. Girko. Spectral Theory of Random Matrices (Russian). Nauka, Moscow, 1988.
[7] Robb J. Muirhead. Aspects of Multivariate Statistical Theory. John Wiley & Sons, New York, 1982.
[8] V.A. Marˇcenko and L.A. Pastur. Distribution of eigenvalues for some sets of random matrices. Math
USSR Sbornik, 1:457–483, 1967.
