Random Matrices

Hsien-Ching Kao

Mathematica Algorithm R&D,  Wolfram Research, Inc.
Champaign, IL  October 20, 2015

Join the Conversation #WolframTechConf

Random matrices

Motivation

Universaility of matrix models in statistics, physics, mathematics, finance, and other engineering fields

Extend the sampling framework in Wolfram Language to handle matrix-valued random variables and their properties

Add support of univariate distributions derived from asymptotic matrix eigenvalues

Applications

Multivariate statistics: regression, PCA, MANOVA, eigenanalysis of covariance, etc.

Matrix normal and matrix t distributions, Wishart distribution, Marchenko–Pastur law

Plot[Evaluate[Table[PDF[MarchenkoPasturDistribution[λ],x],{λ,{0.2,0.5,0.8}}]],{x,0,4},Exclusions->None,Filling->Axis,PlotLabel->Style[" Marchenko–Pastur law",Black,16]]

Physics: nuclear physics, statistical mechanics, quantum chaos, quantum gas, etc.

Gaussian and Circular ensembles, Wigner semicircle and Tracy-Widom laws

GraphicsRow[{Plot[PDF[WignerSemicircleDistribution[1],x],{x,-1,1},Exclusions->None,Filling->Axis,PlotLabel->Style["Wigner semicircle law",Black,16]],Plot[Evaluate[PDF[TracyWidomDistribution[#],x]&/@{1,2,4}],{x,-5,2},Exclusions->None,Filling->Axis,PlotLabel->Style["Tracy-Widom law",Black,16],Epilog->{Text["β=1",{1,0.12}],Text["β=2",{-1,0.39}],Text["β=4",{-1.8,0.55}]}]},ImageSize->700]

Gaussian Ensembles

Gaussian structured matrix distributions

Probability density: P(Z)∝exp(-((Tr[Z.Z])/(2(σ^2))))

GaussianUnitaryMatrixDistribution (GUE): (SuperDagger[x])==x

GaussianUnitaryMatrixDistribution[n]
GaussianUnitaryMatrixDistribution[σ,n]

GaussianOrthogonalMatrixDistribution (GOE): (Transpose[x])==x

GaussianOrthogonalMatrixDistribution[n]
GaussianOrthogonalMatrixDistribution[σ,n]

GaussianSymplecticMatrixDistribution (GSE):  -J.(Transpose[x]).J ==x (self-dual), with J=(({{0, -(Subscript[I, n])}, {(Subscript[I, n]), 0}}))

GaussianSymplecticMatrixDistribution[n]
GaussianSymplecticMatrixDistribution[σ,n]

Matrix property Distribution

An object represents the distribution of lower dimensional properties of matrix-valued random variables

MatrixPropertyDistribution[expr,x\[Distributed]mdist]
MatrixPropertyDistribution[expr,{x1\[Distributed]mdist1, x2\[Distributed]mdist2,...}]

Properties and examples

Wigner semicircle law

n=300;
dist=MatrixPropertyDistribution[Eigenvalues[x],x\[Distributed]GaussianOrthogonalMatrixDistribution[n]];
eigvs=RandomVariate[dist,50];

Show[Histogram[(Join@@eigvs)/Sqrt[2n],{0.025},PDF],Plot[PDF[WignerSemicircleDistribution[1],x],{x,-1,1}]]

Spacing of eigenvalues (Wigner surmise)

dist=MatrixPropertyDistribution[Abs[Subtract@@Eigenvalues[x]],x\[Distributed]GaussianUnitaryMatrixDistribution[2]];gaps=Normalize[RandomVariate[dist,10000],Mean];

Compare sample histogram with Wigner surmise (β=2)

WignerSurmisePDF[x_,β:2]:=2(4 x/Pi)^2 Exp[(-4/Pi) x^2]
Show[Histogram[gaps,{0.1},PDF],Plot[WignerSurmisePDF[x,2],{x,0,2.5}]]

Tracy-Widom law (scaled largest eigenvalue)

n=128;SeedRandom[2];
le=RandomVariate[MatrixPropertyDistribution[(Max[Eigenvalues[x]]-2Sqrt[2.n])(n/4)^(1/6),x\[Distributed]GaussianSymplecticMatrixDistribution[n]],512];

Compare sample histogram with Tracy-Widom law (β=4)

Show[Histogram[le,Automatic,PDF],Plot[PDF[TracyWidomDistribution[4],x ],{x,-5,2},PlotStyle->Thick]]

Circular Ensembles

Uniform distributions on matrix groups

CircularUnitaryMatrixDistribution (CUE): U(n)

CircularUnitaryMatrixDistribution[n]

CircularRealMatrixDistribution (CRE): O(n)

CircularRealMatrixDistribution[n]

CircularQuarternionMatrixDistribution (CQE): USp(2n), x∈U(2n) and  (Transpose[x]).J.x==J

CircularQuaternionMatrixDistribution[n]

Uniform distributions on subsets of unitary group

CircularOrthogonalMatrixDistribution (COE): x∈U(n), (Transpose[x])==x (symmetric)

CircularOrthogonalMatrixDistribution[n]

CircularSymplecticMatrixDistribution (CSE): x∈U(2n), -J.(Transpose[x]).J ==x (self-dual)

CircularSymplecticMatrixDistribution[n]

Properties and examples

Distribution of eigenvalues

dim=100;
args=RandomVariate[MatrixPropertyDistribution[Arg[Eigenvalues[x]],x\[Distributed]CircularUnitaryMatrixDistribution[dim]],100];

Histogram[Join@@args,{-Pi,Pi,2Pi/40},PDF]

Spacing distribution (Wigner surmise)

diffs=Join@@Map[Composition[Differences,Sort],(args/(2Pi))dim];
WignerSurmisePDF[x_,β:2]:=2(4 x/Pi)^2 Exp[(-4/Pi) x^2]
Show[Histogram[diffs,{0.1},PDF],Plot[WignerSurmisePDF[x,2],{x,0,2.5}]]

Boltzmann distribution of Dyson’s Coulomb gas on a circle

Hamiltonian (without kinetic terms): H(Θ)=-(∑_(1≤j<k≤n))log|(e^(i (Subscript[θ, j])))-(e^(i (Subscript[θ, k])))|

hamiltonian[Z_]:=Module[{n=Length[Z]},
Log[n]/2-Sum[Total[Log[Abs[Differences[Z,1,k]]]],{k,n-1}]/n]

Boltzmann distribution: P(Θ)d Θ ∝ exp(-β H(Θ)) dΘ

h𝒟=MatrixPropertyDistribution[hamiltonian[Eigenvalues[x]],x\[Distributed]CircularUnitaryMatrixDistribution[n]];

Compute the sample mean of the Hamiltonian for systems of different size

means=ParallelTable[{n,Mean[RandomVariate[h𝒟,1000]]},{n,10,60,10}]

Plot the sample means and compare them with thermodynamic limit

Show[ListPlot[means,PlotStyle->{Red,PointSize[Medium]},Filling->Axis,PlotRange->{0,0.22},AxesLabel->{"N","⟨H⟩"}],Plot[((1-EulerGamma)/2),{x,0,60}]]

Random rotation

RandomVariate[MatrixPropertyDistribution[r.{0,0,1},r\[Distributed]CircularRealMatrixDistribution[3]]]

Distribution of points over the sphere

Point[RandomVariate[MatrixPropertyDistribution[r.UnitVector[3,3],r\[Distributed]CircularRealMatrixDistribution[3]],10^3]]//Graphics3D

Sample points on (S^3) by randomly rotating a fixed 4D vector

p4d=RandomVariate[MatrixPropertyDistribution[r.UnitVector[4,4],r\[Distributed]CircularRealMatrixDistribution[4]],5 10^4];

Project the points to (S^2)×(S^1) by Hopf map

hopf={x,y,z,w}|->{(ArcTan[x,y]+Pi)/(2Pi),2(x z+y w),2(y z-x w),x^2+y^2-z^2-w^2};

Project the points and bin them by the first coordinate of the projection

proj=hopf@@@p4d;
res=GatherBy[SortBy[proj,First],Quotient[First[#],0.05]&];
res=Drop[res,None,None,1];

Animate[Graphics3D[{PointSize[Small],Point[res[[k]]]},PlotRange->{{-1,1},{-1,1},{-1,1}}],{{k,1,"ϕ"},1,Length[res],1},TrackedSymbols:>k,SaveDefinitions->True,AnimationRunning->False]

Brownian motion in circular ensemble

Construct Brownian motion on CUE by matrices from GUE

mats=RandomVariate[GaussianUnitaryMatrixDistribution[0.1,2],100000];

Generate a Brownian path

initial=RandomVariate[CircularUnitaryMatrixDistribution[2]];
res=FoldList[MatrixExp[I #2].#1&,initial,mats];

Phase of eigenvalues

phases=RandomSample/@Arg[Eigenvalues/@res];

Show[Histogram3D[phases,{-Pi,Pi,0.25},PDF],
Plot3D[(1/(8(Pi^2)))Abs[Exp[I ϕ1]-Exp[I ϕ2]]^2,{ϕ1,-Pi,Pi},{ϕ2,-Pi,Pi},PlotStyle->None,MeshStyle->Thick]]

Matrix-valued Distributions

Wishart and inverse Wishart distributions

Wishart distribution: distribution of empirical covariance matrix

Z=(Transpose[X]).X, X=((Subscript[X, 1]),(Subscript[X, 2]),...), (Subscript[X, k])\[Distributed]N(0,Σ)

with ν degree of freedom

WishartMatrixDistribution[ν,Σ]

Inverse Wishart distribution: inverse of Wishart matrix

InverseWishartMatrixDistribution[ν,Σ]

Properties and examples

Generalization of chi-squared distribution

For any nonzero vector y and Wishart matrix w with scale matrix Σ, y.w.y/(y.Σ.y) is (χ^2) distributed

Σ={{1,4/5},{4/5,1}};
n=10;
y=#/Sqrt[#.Σ.#]&[RandomReal[1,2]];
SeedRandom[1];
data=RandomVariate[MatrixPropertyDistribution[y.w.y,w\[Distributed]WishartMatrixDistribution[n,Σ]],10^4];

DistributionFitTest[data,ChiSquareDistribution[n]]

Show[Histogram[data,Automatic,PDF],Plot[PDF[ChiSquareDistribution[n],x],{x,0,30}]]

Marchenko-Pastur law

eigvsDist[n_,p_]:=MatrixPropertyDistribution[Eigenvalues[x]/n,x\[Distributed]WishartMatrixDistribution[n,IdentityMatrix[p]]]

The limiting distribution of eigenvalues follows MarchenkoPasturDistribution:

data=RandomVariate[eigvsDist[1000,200]];

edist=EstimatedDistribution[data,MarchenkoPasturDistribution[λ,σ]]

Compare the histogram of the eigenvalues with the PDF

Show[Histogram[data,Automatic,PDF],Plot[PDF[edist,x],{x,0,Max[data]+1},PlotStyle->Thick,Exclusions->None]]

Tracy-Widom law for Wishart ensemble

λmax𝒟[n_,p_]:=With[{a=((((Sqrt[n-1/2])+(Sqrt[p-1/2])))^2),b=((Sqrt[n])+(Sqrt[p]))((((1/(Sqrt[n-1/2]))+(1/(Sqrt[p-1/2]))))^(1/3))},MatrixPropertyDistribution[(1/b)(Max[Eigenvalues[𝓂]]-a),𝓂\[Distributed]WishartMatrixDistribution[n,IdentityMatrix[p]]]]

Sample the scaled largest eigenvalue

data=RandomVariate[λmax𝒟[1000,512],250];

Check goodness of fit with TracyWidomDistribution

DistributionFitTest[data,TracyWidomDistribution[1]]

Show[Histogram[data,Automatic,PDF],Plot[PDF[TracyWidomDistribution[1],x],{x,-5,3},PlotStyle->Thick]]

Matrix normal and t distributions

Matrix normal distribution

MatrixNormalDistribution[(Subscript[Σ, row]),(Subscript[Σ, col])]
MatrixNormalDistribution[μ,(Subscript[Σ, row]),(Subscript[Σ, col])]

Matrix t distribution: distribution of MatrixNormalDistribution[Σ,(Subscript[Σ, col])] with Σ sampled from InverseWishartMatrixDistribution[ν+n-1,(Subscript[Σ, row])]

MatrixTDistribution[(Subscript[Σ, row]),(Subscript[Σ, col]),ν]
MatrixTDistribution[μ,(Subscript[Σ, row]),(Subscript[Σ, col]),ν]

Properties and examples

Visualize sample matrices from matrix normal distributions

sigR=Covariance[ARProcess[{-0.95},1][Range[0,100]]];
sigC=KroneckerProduct[IdentityMatrix[100,SparseArray],{{1,1/2},{1/2,3}}];
vals=RandomVariate[MatrixNormalDistribution[sigR,sigC]];

MatrixPlot[vals,PlotTheme->"Minimal"]

sigR=Covariance[ARProcess[{0.95},1][Range[0,100]]];
sigC=KroneckerProduct[IdentityMatrix[100,SparseArray],{{1,1/2},{1/2,3}}];
vals=RandomVariate[MatrixNormalDistribution[sigR,sigC]];

MatrixPlot[vals,PlotTheme->"Minimal"]

Use matrix normal distribution to sample from scalar autoregressive process driven by a multivariate noise signal

sigR=Covariance[ARProcess[{1/2},1][Range[0,10]]];
sigC={{1,1/2},{1/2,3}};
vals=RandomVariate[MatrixNormalDistribution[sigR,sigC],1000];

Construct TemporalData with sampled values and estimate diagonal vector autoregressive process

td=TemporalData[vals,{0,10,1},ValueDimensions->2]

EstimatedProcess[td,ARProcess[{0,0},{a1 IdentityMatrix[2]},{{s11,s12},{s12,s22}}]]

Generalization of t-distribution

For matrix m sampled from matrix t distribution, (Subscript[v, 1]).m.(Subscript[v, 2]) follows Student t distribution

v1=RandomReal[1,2];v2=RandomReal[1,3];

Use MatrixPropertyDistribution to sample values of the expression (Subscript[v, 1]).m.(Subscript[v, 2])

(Subscript[Σ, r])={{1,4/5},{4/5,1}};(Subscript[Σ, c])={{2,1,0},{1,3,1},{0,1,2}};df=5;
data=RandomVariate[MatrixPropertyDistribution[v1.𝓂.v2,𝓂\[Distributed]MatrixTDistribution[(Subscript[Σ, r]),(Subscript[Σ, c]),df]],10^4];

Check agreement with the expected t distribution

t𝒟=StudentTDistribution[0,Sqrt[(v1.(Subscript[Σ, r]).v1)(v2.(Subscript[Σ, c]).v2)/df],df];
DistributionFitTest[data,t𝒟]

Show[Histogram[data,Automatic,PDF],Plot[PDF[t𝒟,x],{x,-4,4},PlotRange->All]]

For matrix m sampled from matrix t distribution, m.v follows multivariate t distribution

v=RandomReal[1,3];

Use MatrixPropertyDistribution to sample values of m.v

(Subscript[Σ, r])={{1,4/5},{4/5,1}};(Subscript[Σ, c])={{2,1,0},{1,3,1},{0,1,2}};df=5;
data=RandomVariate[MatrixPropertyDistribution[𝓂.v,𝓂\[Distributed]MatrixTDistribution[(Subscript[Σ, r]),(Subscript[Σ, c]),df]],10^4];

Verify goodness of fit with the expected distribution

DistributionFitTest[data,MultivariateTDistribution[(Subscript[Σ, r]) (v.(Subscript[Σ, c]).v)/df,df]]

Asymptotic Distributions of Eigenvalues

Marchenko–Pastur distribution

Asymptotic distribution of singular values from Wishart matrices (Z=(X^T).X, X∈M(n,m)  and λ~m/n)

MarchenkoPasturDistribution[λ]
MarchenkoPasturDistribution[λ,σ]

PDF

Plot[Evaluate[Table[PDF[MarchenkoPasturDistribution[λ],x],{λ,{0.2,0.5,0.8}}]],{x,0,4},Exclusions->None,Filling->Axis]

CDF

Plot[Evaluate[Table[CDF[MarchenkoPasturDistribution[λ],x],{λ,{1.0,2.0,4.0}}]],{x,0,8},Filling->Axis,Exclusions->None]

The atomic weight occurs when the Wishart matrix is singular

matrix=Transpose[#].#&[RandomVariate[NormalDistribution[],{500,1000}]];
eigvs=Chop[Eigenvalues[matrix]/500];

Fit MarchenkoPasturDistribution to the eigenvalues

edist=EstimatedDistribution[eigvs,MarchenkoPasturDistribution[λ,1]]

Compare the cumulative histogram of the eigenvalues with the CDF:

Show[Histogram[eigvs,20,CDF],Plot[CDF[edist,x],{x,0,6},PlotStyle->Thick]]

Tracy-Widom distribution

Distribution of scaled largest eigenvalue of random matrices (Gaussian, Wishart ensembles)

PDF

Plot[Evaluate[Table[PDF[TracyWidomDistribution[β],x],{β,{1,2,4}}]],{x,-5,2},Filling->Axis,PlotLegends->{"β = 1","β = 2","β = 4"}]

CDF

Plot[Evaluate[Table[CDF[TracyWidomDistribution[β],x],{β,{1,2,4}}]],{x,-5,2},Filling->Axis,PlotLegends->{"β = 1","β = 2","β = 4"}]

Moments

skewness=Table[Skewness[TracyWidomDistribution[beta]],{beta,{1,2,4}}]

kurtosis=Table[Kurtosis[TracyWidomDistribution[beta]],{beta,{1,2,4}}]

Find their machine-precision approximation

N[skewness]

N[kurtosis]

Compute the Mean of Tracy–Widom distribution to 50-digit precision

N[Mean[TracyWidomDistribution[2]],50]

Tail of CDF

(Subscript[F, β_])[x_]:=CDF[TracyWidomDistribution[β],x];
logCDFapprox1[x_]=-((11Log[2])/48)+(1/24)-(1/2)Log[Glaisher]-(((Abs[x])^3)/24)-(((Abs[x])^(3/2))/(3(Sqrt[2])))-(1/16)Log[Abs[x]];
logCDFapprox2[x_]=((Log[2])/24)+(1/12)-Log[Glaisher]-(((Abs[x])^3)/12)-(1/8)Log[Abs[x]];
logCDFapprox4[x_]=-((73 Log[2])/96)+(1/24)-(1/2)Log[Glaisher]-(((Abs[x])^3)/(6(Sqrt[2])))+(((2^(1/4))((Abs[x])^(3/2)))/3)-(1/16)Log[Abs[x]];

Show[LogPlot[{(Subscript[F, 1])[x],(Subscript[F, 2])[x],(Subscript[F, 4])[x]},{x,-12,-5},Filling->Axis,PlotLegends->"Expressions"],Plot[{logCDFapprox1[x],logCDFapprox2[x],logCDFapprox4[x]},{x,-12,-5},PlotStyle->Directive[Black,Dotted]],PlotLabel->"Tracy-Widom CDF in Log scale"]

Plot[{Exp[logCDFapprox1[x]]/(Subscript[F, 1])[x],Exp[logCDFapprox2[x]]/(Subscript[F, 2])[x],Exp[logCDFapprox4[x]]/(Subscript[F, 4])[x]},{x,-12,-5},PlotStyle->Automatic,PlotLegends->{"(Subscript[F, 1])(x)","(Subscript[F, 2])(x)","(Subscript[F, 4])(x)"},PlotLabel->"Relative error"]

Approximate Tracy–Widom distribution by GammaDistribution

gdist=GammaDistribution[k,s,1,a];
PDF[gdist,x]

Match GammaDistribution with Tracy–Widom distribution (β=2) with the leading moments:

moments=Through[{Mean,Variance,Skewness}[gdist]];
nmoments=N[Through[{Mean,Variance,Skewness}[TracyWidomDistribution[2]]]];
Quiet[sol=FindRoot[Thread[moments==nmoments],{{k,1},{s,1},{a,1}}]]

Compare the PDFs between -5≤x≤2:

Plot[{Evaluate[PDF[gdist/.sol,x]],PDF[TracyWidomDistribution[2],x]},{x,-5,2},PlotLegends->{"Gamma distribution","Tracy-Widom distribution"}]

Compare the CDFs between a+1≤x≤-5 in log scale:

LogPlot[{Evaluate[CDF[gdist/.sol,x]],CDF[TracyWidomDistribution[2],x]},{x,a+1/.sol,-5},PlotLegends->{"Gamma distribution","Tracy-Widom distribution"}]

Distribution of longest increasing subsequence

losLength[data_]:=Length[LongestOrderedSequence[data]]

Find lengths of the longest increasing subsequence in random permutations of list {1,2,…,n}:

n=50000;list=Range[n];
lengths=ParallelTable[losLength[RandomSample[list]],{2000}];

Compare the scaled length of the longest increasing subsequence with Tracy–Widom distribution of β=2:

scaledLength=(1/(n^(1/6)))(lengths-2.0(Sqrt[n]));

Show[SmoothHistogram[scaledLength,PlotLegends->{"Scaled length"}],Plot[PDF[TracyWidomDistribution[2],x],{x,-4,2},PlotStyle->,PlotLegends->{"Tracy-Widom distribution"}]]

Final Remark

Sampling of matrix-valued random variables

Gaussian ensembles

Circular ensembles

Matrix-valued distributions (Wishart, Normal, t distribution)

Sampling of properties of matrix-valued random variables

MatrixPropertyDistribution

New matrix related univariate distributions

Marchenko-Pastur law

Tracy-Widom law

Thank you for your attention!