Numerical%20Computations%20and%20Random%20Matrix%20Theory

1
Numerical Computations and Random Matrix Theory

• Alan Edelman
• MIT Dept of Mathematics,
• Computer Science AI Laboratories
• Friday February 25, 2005

2
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')
3
Tidbit of interest to Matrix Computations
Audience

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

4
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

5
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
6
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
7
Tidbit of interest to Matrix Computations
Audienceand pure mathematicians!

• The most analytical random matrices seen from on
  high

8
Same structure everywhere!
Orthog Matrix MATLAB (Arandn(n)
Brandn(n))
Hermite Sym Eig eig(AA)
Laguerre SVD eig(AA)
Jacobi GSVD gsvd(A,B)
Fourier Eig U,Rqr(AiB)
9
Same structure everywhere!
Orthog Matrix Weight Stats
Graph Theory SymSpace
Hermite Sym Eig exp(-x2) Normal Complete Graph non-compact A,AI,AII
Laguerre SVD xae-x Chi-squared Bipartite Graph non-compact AIII,BDI,CII
Jacobi GSVD (1-x)a x (1x)ß Beta Regular Graph compact A, AI, AII, C, D, CI, D, DIII
Fourier Eig ei? compact AIII, BDI, CDI
10
Tidbit of interest to Matrix Computations
Audienceand combinatorists!

• The longest increasing subsequence

11
Longest Increasing Subsequence(n4)
Green 4 Yellow 3 Red 2 Purple 1
1 2 3 4 2 1 3 4 3 1 2 4 4 1 2 3
1 2 4 3 2 1 4 3 3 1 4 2 4 1 3 2
1 3 2 4 2 3 1 4 3 2 1 4 4 2 1 3
1 3 4 2 2 3 4 1 3 2 4 1 4 2 3 1
1 4 2 3 2 4 1 3 3 4 1 2 4 3 1 2
1 4 3 2 2 4 3 1 3 4 2 1 4 3 2 1
12
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

13
Haar or not Haar?
14
Tidbit!

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

15
sym matrix to tridiagonal form
G ?6
?6 G ?5
?5 G ?4
?4 G ?3
?3 G ?2
?2 G ?1
?1 G
Same eigenvalue distribution as AA O(n)
storage !! O(n) compute
16
General beta
G ?6?
?6? G ?5?
?5? G ?4?
?4? G ?3?
?3? G ?2?
?2? G ??
?? G
beta 1 reals 2 complexes 4 quaternions
Bidiagonal Version corresponds To Wishart
matrices of Statistics
17
Largest Eigenvalue of Hermite
18
Painlevé Equations
19
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)

20
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.)

21
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

22
Spacings of eigs of AA
23
Riemann Zeta Zeros
24
Stochastic Operator
25
Everyones Favorite Tridiagonal
-2 1
1 -2 1

1
1 -2





26
Everyones Favorite Tridiagonal
-2 1
1 -2 1

1
1 -2
G
G


G
1 (ßn)1/2







27
Stochastic Operator Limit



28
Tidbit

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

29
Free Probability vs Classical Probability
30
Random Matrix Calculator
31
How to use calculator
32
Steps 1 and 2
33
Steps 3 and 4
34
Steps 5 and 6
35
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

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

• Linear Algebra Randomness !!!

42
Mops (Dumitriu etc.) Symbolic
43
Symbolic MOPS applications
ß3 hist(eig(S))
44
(No Transcript)
45
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

46
Largest Eigenvalue Plots
47
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)

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

• The longest increasing subsequence

51
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???

52
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!

?
53
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 ?
54
Random cond numbers, n??
Distribution of ?/n
Experiment with n200
55
Finite n

• n10
  n25
• n50
  n100

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

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

Complex n x n, n?8
57
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