Numerical Computations and Random Matrix Theory

1
Numerical Computations and Random Matrix Theory

• Alan Edelman
• MIT Dept of Mathematics,
• Computer Science AI Laboratories
• Tuesday May 24, 2005

2
Applications

• 65 Signal Processing People Signed up for
  tutorial
• \ eig() Stochastic Differential Equations
  Stochastic Eigenanalysis

3
Wigners Semi-Circle

• The classical most famous rand eig theorem
• Let S random symmetric Gaussian
• MATLAB Arandn(n) S(AA)/2
• S known as the Gaussian Orthogonal Ensemble
• Normalized eigenvalue histogram is a semi-circle
• Precise statements require n?? etc.

n20 s30000 d.05 matrix size, samples,
sample dist e gather up
eigenvalues im1 imaginary(1) or
real(0) for i1s, arandn(n)imsqrt(-1)randn(
n)a(aa')/(2sqrt(2n(im1))) veig(a)'
ee v end hold off m xhist(e,-1.5d1.5)
bar(x,mpi/(2dns)) axis('square') axis(-1.5
1.5 -1 2) hold on t-1.011
plot(t,sqrt(1-t.2),'r')
4
Matrix Calculus????

• Condition Numbers and Jacobians of Matrix
  Functions and Factorizations or
• What is matrix calculus??

5
Matrix Functions and Factorizations

• e.g. f(A)A2 or L,Ulu(A) or Q,Rqr(A)
• U, R n(n1)/2 parameters
• L, Q n(n-1)/2 parameters
• Q globally (Householder)
• Q locally (tangent space Qantisym )
• The Jacobian or df or linearization is n2
  x n2
• fS?S2 (sym) df is n(n1)/2 x n(n1)/2
• fQ?Q2 (orth) df is n(n-1)/2 x n(n-1)/2

6
Condition number of a matrix function or
factorization
Jacobian Det J ?si(df)det(df) Example 1
f(A)A2 df(A) kron(I,A)kron(AT ,I) Example 2
f(A)A-1 df(A)-kron(A-T,A-1) df(A)A-12
? ?A A-1
7
Matrix Factorization Jacobians
General
ALU AU?VT AX?X-1
AQR AQS (polar)
? uiin-i
? riim-i
? (?i2- ?j2)
? (?i?j)
? (?i-?j)2
Sym
Orthogonal
?sin(?i ?j)sin (?i- ?j)
Tridiagonal
TQ?QT
? (ti1,i)/ ?qi
8
Tidbit of interest to Matrix Computations
Audienceand pure mathematicians!

• The most analytical random matrices seen from on
  high

9
Same structure everywhere!
Orthog Matrix MATLAB (Arandn(n)
Brandn(n))
10
Same structure everywhere!
Orthog Matrix Weight Stats
Graph Theory SymSpace
11
Tidbit of interest to Matrix Computations
Audienceand combinatorists!

• The longest increasing subsequence

12
Longest Increasing Subsequence(n4)
Green 4 Yellow 3 Red 2 Purple 1
13
Random Matrix Result

• Permutations on 1..n with longest increasing
  subsequence k is
• E ( tr(Qk)2n) . The 2nth moment of the
  absolute trace of random kxk orthogonal matrices
• Longest increasing subsequence is the parallel
  complexity of an upper triangular solve with
  sparsity given by
• Uij(p) ?0 if p(i)p(j) and ij

14
Haar or not Haar?
15

• Random Tridiagonalization leads to eigenvalues of
  billion by billion matrix!

16
Largest Eigenvalue of Hermite
17
Painlevé Equations
18
MATLAB

• beta1 n1e9 opts.disp0opts.issym1
• alpha10 kround(alphan(1/3)) cutoff
  parameters
• dsqrt(chi2rnd( beta(n-1(n-k-1))))'
• Hspdiags( d,1,k,k)spdiags(randn(k,1),0,k,k)
• H(HH')/sqrt(4nbeta)
• eigs(H,1,1,opts)

19
Tricks to get O(n9) speedup

• Sparse matrix storage (Only O(n) storage is used)
  
• Tridiagonal Ensemble Formulas (Any beta is
  available due to the tridiagonal ensemble)
• The Lanczos Algorithm for Eigenvalue Computation
  ( This allows the computation of the extreme
  eigenvalue faster than typical general purpose
  eigensolvers.)
• The shift-and-invert accelerator to Lanczos and
  Arnoldi (Since we know the eigenvalues are near
  1, we can accelerate the convergence of the
  largest eigenvalue)
• The ARPACK software package as made available
  seamlessly in MATLAB (The Arnoldi package
  contains state of the art data structures and
  numerical choices.)
• The observation that if k 10n1/3 , then the
  largest eigenvalue is determined numerically by
  the top k k segment of n. (This is an
  interesting mathematical statement related to the
  decay of the Airy function.)

20
Tidbit of interest to Matrix Computations
AudienceStochastic Eigenequations

• Continuous vs Discrete
• Diff Eqns Matrix Comps Cont Eig Matrix
  Eigs
• Add probability
• Stochastic Differential Equations Stochastic
  Eigenequations
• Finite Random Matrix Theory

21
Spacings of eigs of AA
22
Riemann Zeta Zeros
23
Stochastic Operator
24
Everyones Favorite Tridiagonal





25
Everyones Favorite Tridiagonal
1 (ßn)1/2







26
Stochastic Operator Limit



27
Tidbit

• eig(AB) eig(A) eig(B) ?????

28
Free Probability vs Classical Probability
29
Random Matrix Calculator
30
How to use calculator
31
Steps 1 and 2
32
Steps 3 and 4
33
Steps 5 and 6
34
Multivariate Orthogonal PolynomialsHypergeometr
ics of Matrix Argument

• The important special functions of the 21st
  century
• Begin with w(x) on I
• ? p?(x)p?(x) ?(x)ß ?i w(xi)dxi d??
• Jack Polynomials orthogonal for w1 on the unit
  circle. Analogs of xm

35
Multivariate Hypergeometric Functions
36
Multivariate Hypergeometric Functions
37
Plamens clever idea
38
Smallest eigenvalue statistics
Arandn(m,n) hist(min(svd(A).2))
39
Symbolic MOPS applications
Arandn(n) S(AA)/2 trace(S4)
det(S3)
40
Summary

• Linear Algebra Randomness !!!

41
Mops (Dumitriu etc.) Symbolic
42
Symbolic MOPS applications
ß3 hist(eig(S))
43
(No Transcript)
44
Spacings

• Take a large collection of consecutive
  zeros/eigenvalues.
• Normalize so that average spacing 1.
• Spacing Function Histogram of consecutive
  differences (the (k1)st the kth)
• Pairwise Correlation Function Histogram of all
  possible differences (the kth the jth)
• Conjecture These functions are the same for
  random matrices and Riemann zeta

45
Largest Eigenvalue Plots
46
MATLAB

• beta1 n1e9 opts.disp0opts.issym1
• alpha10 kround(alphan(1/3)) cutoff
  parameters
• dsqrt(chi2rnd( beta(n-1(n-k-1))))'
• Hspdiags( d,1,k,k)spdiags(randn(k,1),0,k,k)
• H(HH')/sqrt(4nbeta)
• eigs(H,1,1,opts)

47
Open Problems
The distribution for general beta Seems to be
governed by a convection-diffusion equation
48
Random matrix tools!
49
Tidbit of interest to Matrix Computations
Audienceand combinatorists!

• The longest increasing subsequence

50
Numerical Analysis Condition Numbers

• ?(A) condition number of A
• If AU?V is the svd, then ?(A) ?max/?min .
• Alternatively, ?(A) ?? max (AA)/?? min (AA)
• One number that measures digits lost in finite
  precision and general matrix badness
• Smallgood
• Largebad
• The condition of a random matrix???

51
Von Neumann co.

• Solve Axb via x (AA) -1A b
• M ?A-1
• Matrix Residual AM-I2
• AM-I2lt 200?2 n ?
• How should we estimate ??
• Assume, as a model, that the elements of A are
  independent standard normals!

?
52
Von Neumann co. estimates (1947-1951)

• For a random matrix of order n the expectation
  value has been shown to be about n
• Goldstine, von Neumann
• we choose two different values of ?, namely n
  and ?10n
• Bargmann, Montgomery, vN
• With a probability 1 ? lt 10n
• Goldstine, von Neumann

X ?
53
Random cond numbers, n??
Distribution of ?/n
Experiment with n200
54
Finite n

• n10
  n25
• n50
  n100

55
Condition Number Distributions
Real n x n, n?8

• Generalizations
• ß 1real, 2complex
• finite matrices
• rectangular m x n

Complex n x n, n?8
56
Condition Number Distributions
Square, n?8 P(?/nß gt x) (2ßß-1/G(ß))/xß
(All Betas!!) General Formula P(?gtx) Cµ/x
ß(n-m1), where µ ß(n-m1)/2th moment of the
largest eigenvalue of Wm-1,n1 (ß) and C
is a known geometrical constant. Density for the
largest eig of W is known in terms of
1F1((ß/2)(n1), ((ß/2)(nm-1) -(x/2)Im-1) from
which µ is available Tracy-Widom law applies
probably all beta for large m,n. Johnstone shows
at least beta1,2.
Real n x n, n?8
Complex n x n, n?8