18.338J/16.394J: The Mathematics of Infinite Random Matrices
Tridiagonal Matrices, Orthogonal Polynomials and the Classical
Random Matrix Ensembles
Brian Sutton
Handout #5, Thursday, September 23, 2004
In class, we saw the connection between the so-called Hermite matrix and the semi-circular law. There
is actually a deeper story that connects the classical random matrix ensembles to the classical orthogonal
polynomials studied in classical texts such as [1] and more recent monographs such as [2]. We illuminate
part of this story here. The website www.mathworld.com is an excellent reference for these polynomials and
will prove handy when completing the exercises.
In any computational explorations, see if you can spot the interesting feature in the eigenvectors (either
the first or last row/column) of the corresponding tridiagonal matrix.
1 Roots of orthogonal polynomials
Any weight function w(x) on an interval [a,b] (possibly infinite) defines a system of orthogonal polynomials
pin, n = 0,1,2,..., satisfying
integraldisplay b
a
pimpinw(x)dx = δmn. (1)
The polynomials can be generated easily, because they satisfy a three-term recurrence
pi0(x) =
parenleftBigintegraltextb
a w
parenrightBig?1/2
, (2)
xpin(x) = bnpin?1(x) + an+1pin(x) + bn+1pin+1(x). (3)
(Note that b0 is taken to be zero.)
Perhaps surprisingly, the three-term recurrence can be used to find the roots of pin as well. Simply form
the symmetric tridiagonal matrix
Tn =
?
??
??
??
??
a1 b1
b1 a2 b2
b2 a3 b3
... ... ...
bn?2 an?1 bn?1
bn?1 an
?
??
??
??
??
.
The roots of pin are the eigenvalues of Tn!
Remark 1. The symmetry of T follows from two choices which we made. First, in (1), pin was normalized to
have unit length under the inner product 〈f,g〉 = integraltextba fgw. Second, (3) was arranged in the form xpin = RHS.
Not all authors follow these conventions.
To see that the eigenvalues are the roots, consider vectors of the form
v = (pi0(x),pi1(x),...,pin?1(x))T.
The recurrence (3) implies that (Tv)i = xpii?1(x) = xvi for i = 1,...,n ? 1, regardless of the value of x.
Now, if x is a root of pin(x), then (3) specializes to
xpin?1(x) = bn?1pin?2(x) + anpin?1x
i.e.
xvn = bn?1vn?1 + anvn = (Tv)n.
Therefore, every root of pin(x) is an eigenvalue of Tn.
2 Hermite polynomials
Let w(x) = e?x2 on R. According to mathworld.wolfram.com, the system of polynomials ?Hn defined by
integraldisplay ∞
?∞
?Hm ?Hnw(x)dx = δmn2nn!√pi
satisfies the recurrence
?H0(x) = 1 (4)
?Hn+1(x) = 2x ?Hn(x)?2n ?Hn?1(x). (5)
Let the orthonormal Hermite polynomials be defined by Hn(x) = 12n/2√n!pi1/4 ?Hn(x), so that
integraldisplay ∞
?∞
HmHnwdx = δmn,
and substitute 2n/2(n!)1/2pi1/4Hn(x) for ?Hn(x) in (5) to find
H0(x) = pi1/4
xHn(x) = 1√2√nHn?1(x) + 1√2√n + 1Hn+1(x).
These equations are exactly of the form (1)-(3), so the eigenvalues of
1√
2
?
??
??
??
??
0 √1√
1 0 √2
0 √2 0 √3
... ... ...
√n?2 0 √n?1
√n?1 0
?
??
??
??
??
are the roots of Hn.
Exericse: Write a MATLAB function that generates this tridiagonal matrix. Use the eig command to
compute its eigenvalues for any choice of n . Histogram the roots of Hn for n = 50,100,250,300,500
using the histn function. Compare it to the semi-circular law. Does it make a difference whether the the
generalized ensemble is real or complex?
3 Laguerre polynomials
Let wα(x) = xαe?x on [0,∞).
Exercise: Verify that the Laguerre polynomials defined by
integraldisplay ∞
0
LαmLαnwdx = δmn
satisfy the recurrence
Lα0(x) = 1radicalbigΓ(α + 1)
xLαn(x) = ?√n√α+ nLαn?1(x) + (α + 1 + 2n)Lαn(x)?√n + 1√α + n + 1Lαn+1(x).
Therefore, the roots of Lαn are given by the eigenvalues of
?
??
??
??
α+ 1 ?√1√α + 1
?√1√α + 1 α + 3 ?√2√α+ 2
0 ?√2√α + 2 α + 5 ?√3√α + 3
... ... ...
?√n?1√α + n?1 α?1 + 2n
?
??
??
??.
Exericse: Write a MATLAB function that generates this tridiagonal matrix. Use the eig command to
compute its eigenvalues for any choice of n and α. Histogram the roots of Ln for different values of n and
m using the histn function. Compare it to the Marˇcenko-Pastur distribution for the generalized Laguerre
Ensemble. Does it make a difference whether the the generalized ensemble is real or complex? Note:
α = β (n?m+1)/2?1 where β = 1 when the elements are real and β = 2 when the elements are complex.
The parameters m and n are associated with the rows and columns of the Gaussian matrix used to
construct the Wishart matrix. Recall that m/n = c.
4 Jacobi polynomials
Jacobi polynomials are orthogonal with respect to the weight w(x) = (1?x)α(1 + x)β on the interval
[?1,1].
Exercise: Verify that the tridiagonal matrix describing the recurrence is produced by the following MATLAB
code.
a39
a38
a36
a37
function t = trijacobi(n,a,b)
%TRIJACOBI Returns the tridiagonal matrix of polynomial recurrence coefficients
% TRIJACOBI(N,A,B) returns a tridiagonal matrix whose entries correspond
% to the coefficients describing the recurrence between
% the Jacobi Polynomials.
%
% The eigenvalues of this matrix correspond to the roots of the N-th
% Jacobi polynomial with parameters A and B.
%
% N, A and B are integers.
%
%
% References:
% [1] Alan Edelman, Handout 6: Tridiagonal Matrices, Orthogonal
% Polynomials ans the Classical Random Matrix
% Ensembles, Fall 2004,
% Course Notes 18.338.
% [2] Alan Edelman, Random Matrix Eigenvalues.
% [3] Gabor Szego, Orthogonal Polynomials, American Mathematical
% Society, Providence, 1975. 4th Edition.
% [4] Eric Weisstein, Jacobi Polynomial." From MathWorld--
% A Wolfram Web Resource.
% http://mathworld.wolfram.com/JacobiPolynomial.htm
%
% Brian Sutton, Sept. 2004.
% $Revision: 1.0 $ $Date: 2004/09/23 11:35:18 $
i = (1:n-1)’;
d1 = 2*sqrt(i.*(a+i).*(b+i).*(a+b+i)./(-1+a+b+2*i)./(1+a+b+2*i))./(a+ ...
b+2*i);
d0 = -(a^2-b^2)./(a+b+2*(0:n-1)’)./(2+a+b+2*(0:n-1)’);
t = spdiags([[d1;0] d0 [0;d1]],[-1 0 1],n,n);
Code 1
Exericse: Use the eig command to compute its eigenvalues for any choice of n. Histogram the roots of Tn
for different values of n, a and b using the histn function. Compare it to the MANOVA density for the
MANOVA matrix. Hint: You might have to shift the spectrum for it to line up exactly. Does it make a
difference whether the the generalized ensemble is real or complex? Note: a = β (n1 ?m + 1)/2?1 and
b = β (n2 ?m + 1)/2?1 where β = 1 when the elements are real and β = 2 when the elements are
complex. The parameters m, n1 and n2 are associated with the rows and columns of the Gaussian matrices
used to construct the MANOVA matrix. Recall that m/n1 = c1 and m/n2 = c2.
References
[1] G′abor Szeg¨o. Orthogonal Polynomials. American Mathematical Society, Providence, 1975. 4th edition.
[2] Percy Deift. Orthogonal Polynomials and Random Matrices: A Riemann-Hilbert Approach. American
Mathematical Society, Providence, 1998.
