Numerical Computations and Random Matrix Theory Alan Edelman

  • Slides: 31
Download presentation
Numerical Computations and Random Matrix Theory Alan Edelman MIT: Dept of Mathematics, Computer Science

Numerical Computations and Random Matrix Theory Alan Edelman MIT: Dept of Mathematics, Computer Science AI Laboratories Tuesday May 24, 2005 9/6/2021 1

Applications v 65 Signal Processing People Signed up for tutorial v “”: eig() :

Applications v 65 Signal Processing People Signed up for tutorial v “”: eig() : : Stochastic Differential Equations : Stochastic Eigenanalysis 9/6/2021 2

Matrix Calculus? ? v Condition Numbers and Jacobians of Matrix Functions and Factorizations or

Matrix Calculus? ? v Condition Numbers and Jacobians of Matrix Functions and Factorizations or v What is matrix calculus? ? 9/6/2021 4

Matrix Functions and Factorizations v e. g. f(A)=A 2 or [L, U]=lu(A) or [Q,

Matrix Functions and Factorizations v e. g. f(A)=A 2 or [L, U]=lu(A) or [Q, R]=qr(A) v U, R: n(n+1)/2 parameters v L, Q: n(n-1)/2 parameters v. Q globally (Householder) v Q locally (tangent space = Q*antisym ) v The Jacobian or “df” or “linearization” is n 2 x n 2 v f: S S 2 (sym) df is n(n+1)/2 x n(n+1)/2 v f: Q Q 2 (orth) df is n(n-1)/2 x n(n-1)/2 9/6/2021 5

Condition number of a matrix function or factorization Jacobian Det = J = ∏σi(df)=det(df)

Condition number of a matrix function or factorization Jacobian Det = J = ∏σi(df)=det(df) Example 1: f(A)=A 2 df(A) =kron(I, A)+kron(AT , I) 9/6/2021 Example 2: f(A)=A-1 df(A)=-kron(A-T, A-1) ||df(A)||=||A-1||2 κ=||A|| ||A-1|| 6

Matrix Factorization Jacobians General uiin-i riim-i A=LU A=QR A=U VT ( i 2 -

Matrix Factorization Jacobians General uiin-i riim-i A=LU A=QR A=U VT ( i 2 - j 2) A=QS (polar) ( i+ j) A=X X-1 ( i- j)2 Sym S=Q QT S=LLT S=LDLT Orthogonal ( i- j) 2 n liin+1 -i Q=U C S VT sin( + )sin ( - ) i j S -C din-i [ ] Tridiagonal T=Q QT (ti+1, i)/ qi 9/6/2021 7

Same structure everywhere! Orthog Matrix MATLAB (A=randn(n) B=randn(n)) Hermite Sym Eig eig(A+A’) Laguerre SVD

Same structure everywhere! Orthog Matrix MATLAB (A=randn(n) B=randn(n)) Hermite Sym Eig eig(A+A’) Laguerre SVD eig(A*A’) Jacobi GSVD gsvd(A, B) Fourier Eig [U, R]=qr(A+i*B) 9/6/2021 9

Same structure everywhere! Orthog Matrix Weight Stats Hermite Sym Eig exp(-x 2) Normal Laguerre

Same structure everywhere! Orthog Matrix Weight Stats Hermite Sym Eig exp(-x 2) Normal Laguerre SVD Jacobi GSVD Fourier Eig 9/6/2021 xαe-x Chisquared (1 -x)α x Beta β (1+x) eiθ Graph Theory Sym. Space Complete Graph Bipartite Graph noncompact A, AII noncompact AIII, BDI, CII compact Regular Graph A, AII, C, D, CI, D, DIII compact AIII, BDI, 10 CDI

Haar or not Haar? 9/6/2021 14

Haar or not Haar? 9/6/2021 14

v Random Tridiagonalization leads to eigenvalues of billion by billion matrix! 9/6/2021 15

v Random Tridiagonalization leads to eigenvalues of billion by billion matrix! 9/6/2021 15

Largest Eigenvalue of Hermite 9/6/2021 16

Largest Eigenvalue of Hermite 9/6/2021 16

MATLAB beta=1; n=1 e 9; opts. disp=0; opts. issym=1; alpha=10; k=round(alpha*n^(1/3)); % cutoff parameters

MATLAB beta=1; n=1 e 9; opts. disp=0; opts. issym=1; alpha=10; k=round(alpha*n^(1/3)); % cutoff parameters d=sqrt(chi 2 rnd( beta*(n: -1: (n-k-1))))'; H=spdiags( d, 1, k, k)+spdiags(randn(k, 1), 0, k, k); H=(H+H')/sqrt(4*n*beta); eigs(H, 1, 1, opts) 9/6/2021 18

Tricks to get O(n 9) speedup • Sparse matrix storage (Only O(n) storage is

Tricks to get O(n 9) 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 = 10 n 1/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. ) 9/6/2021 19

Spacings of eigs of A+A’ 9/6/2021 21

Spacings of eigs of A+A’ 9/6/2021 21

Riemann Zeta Zeros 9/6/2021 22

Riemann Zeta Zeros 9/6/2021 22

Stochastic Operator 9/6/2021 23

Stochastic Operator 9/6/2021 23

Everyone’s Favorite Tridiagonal 1 n 2 -2 1 1 -2 1 … … …

Everyone’s Favorite Tridiagonal 1 n 2 -2 1 1 -2 1 … … … 1 1 -2 d 2 dx 2 9/6/2021 24

Everyone’s Favorite Tridiagonal 1 n 2 -2 1 1 -2 1 … … …

Everyone’s Favorite Tridiagonal 1 n 2 -2 1 1 -2 1 … … … 1 1 -2 G d 2 dx 2 9/6/2021 G 1 +(βn)1/2 G + d. W β 1/2 25

Tidbit v eig(A+B) = eig(A) + eig(B) ? ? ? 9/6/2021 27

Tidbit v eig(A+B) = eig(A) + eig(B) ? ? ? 9/6/2021 27

Free Probability vs Classical Probability 9/6/2021 28

Free Probability vs Classical Probability 9/6/2021 28

Random Matrix Calculator 9/6/2021 29

Random Matrix Calculator 9/6/2021 29

How to use calculator 9/6/2021 30

How to use calculator 9/6/2021 30

Steps 1 and 2 9/6/2021 31

Steps 1 and 2 9/6/2021 31

Steps 3 and 4 9/6/2021 32

Steps 3 and 4 9/6/2021 32

Steps 5 and 6 9/6/2021 33

Steps 5 and 6 9/6/2021 33

Multivariate Orthogonal Polynomials & Hypergeometrics of Matrix Argument v The important special functions of

Multivariate Orthogonal Polynomials & Hypergeometrics of Matrix Argument v The important special functions of the 21 st century v Begin with w(x) on I v∫ pκ(x)pλ(x) Δ(x)β ∏i w(xi)dxi = δκλ v. Jack Polynomials orthogonal for w=1 on the unit circle. Analogs of xm 9/6/2021 34

Multivariate Hypergeometric Functions 9/6/2021 35

Multivariate Hypergeometric Functions 9/6/2021 35

Multivariate Hypergeometric Functions 9/6/2021 36

Multivariate Hypergeometric Functions 9/6/2021 36

Plamen’s clever idea 9/6/2021 37

Plamen’s clever idea 9/6/2021 37

Summary v Linear 9/6/2021 Algebra + Randomness !!! 40

Summary v Linear 9/6/2021 Algebra + Randomness !!! 40