18.325 – Random Matrices Numerical Methods for Random Matrices Per-Olof Persson ([email protected]) December 19, 2002� � 1 Largest Eigenvalue Distributions In this section, the distributions of the largest eigenvalue of matrices in the β-ensembles are studied. Histograms are created first by simulation, then by solving the Painlev´e II nonlinear differential equation. 1.1 Simulation The Gaussian Unitary Ensemble (GUE) is defined as the Hermitian n × n matrices A, where the diagonal elements xjj and the upper triangular elements xjk = ujk + ivjk are independent Gaussians with zero-mean, and Var(xjj) = 1, 1 ≤ j ≤ n, (1)1Var(ujk) = Var(vjk) = 2, 1 ≤ j < k ≤ n. Since a sum of Gaussians is a ne w Gaussian, an instance of these matrices can be created conveniently in MATLAB: A=randn(n)+i*randn(n); A=(A+A’)/2; The largest eigenvalue of this matrix is about 2√n. To get a distribution that converges as n → ∞, the shifted and scaled largest eigenvalue λ�ismax calculated as λ�= nmax 1 6λmax − 2√n � . (2) It is now straight-forward to compute the distribution for λ�by simula-max tion: for ii=1:trials A=randn(n)+i*randn(n); A=(A+A’)/2; lmax=max(eig(A)); lmaxscaled=n^(1/6)*(lmax-2*sqrt(n)); % Store lmax end % Create and plot histogram The problem with this technique is that the computational requirements and the memory requirements grow fast with n, which should be as large as possible in order to be a good approximation of infinity. Just storing the matrix A requires n2 double-precision numbers, so on most computers 13 today n has to be less than 104. Furthermore, computing all the eigenvalues of a full Hermitian matrix requires a computing time proportional to n . This means that it will take many days to create a smooth histogram by simulation, even for relatively small values of n. To improve upon this situation, another matrix can be studied that has the s ame eigenvalue distribution as A above. In [2], it was shown that this is true for the following symmetric matrix when β = 2: ⎞⎛ N(0, 2) χ(n−1)β χ(n−1)β N(0, 2) χ(n−2)β1 Hβ ∼ √2 ⎜⎜⎜⎜⎜⎝ ⎟⎟⎟⎟⎟⎠ . (3) . . . . . .. . . χ2β N(0, 2) χβ χβ N(0, 2) Here, N(0, 2) is a zero-mean Gaussian with variance 2, and χd is the square-root of a χ2 distributed number with d degrees of freedom. Note that the matrix is symmetric, so the subdiagonal and the superdiagonal are always equal. This matrix has a tridiagonal sparsity structure, and only 2n double-precision numbers are required to store an instance of it. The time for computing the largest eigenvalue is proportional to n, either using Krylov subspace based methods or the method of bisection [7]. In MATLAB, the built-in function eigs can be used, although that re-quires dealing with the sparse matrix structure. There is also a large amount of overhead in this function, which results in a relatively poor performance. Instead, the function maxeig is used below to compute the eigenvalues. This is not a built-in function in MATLAB, but it can be downloaded from http://www.mit.edu/~persson/mltrid. It is based on the method of bi-section, and requires just two ordinary MATLAB vectors as input, corre-1sponding to the diagonal and the subdiagonal. It also turns out that only the first 10n components of the eigenvector 3 corresponding to the larges t eigenvalue are significantly greater than zero. Therefore, the upper-left ncutoff by ncutoff submatrix has the same largest eigenvalue (or at least very close), where ncutoff ≈ 10n 1 3. (4) Matrices of size n = 1 1012 can then easily be used since the computations can 105. Also, for these large values of be done on a matrix of size only 10n 3= n the approximation χ2 n ≈ n is accurate. 2A histogram of the distribution for n = 109 can now be created using the code below. n=1e9; nrep=1e4; beta=2; cutoff=round(10*n^(1/3)); d1=sqrt(n-1:-1:n+1-cutoff)’/2/sqrt(n); ls=zeros(1,nrep); for ii=1:nrep d0=randn(cutoff,1)/sqrt(n*beta); ls(ii)=maxeig(d0,d1); end ls=(ls-1)*n^(2/3)*2; histdistr(ls,-7:0.2:3) where the function histdistr below is used to histogram the data. It assumes that the histogram boxes are equidistant. function [xmid,H]=histdistr(ls,x) dx=x(2)-x(1); H=histc(ls,x); H=H(1:end-1); H=H/sum(H)/dx; xmid=(x(1:end-1)+x(2:end))/2; bar(xmid,H) grid on The resulting distribution is shown in Figure 1, together with distribu-tions for β = 1 and β = 4. The plots also contain solid curves representing the “true solutions” (see next section). 1.2 Painlev´e II Instead of using simulation to plot the distributions of the largest eigen-values, it can be computed from the solution of the Painlev´e II nonlinear differential equation [6]: 3 q�� = sq + 2q (5) 3� � � � � −6 −4 −2 0 2 400.10.20.30.40.50.6Normalized and Scaled Largest EigenvalueProbabilityβ=1−6 −4 −2 0 2 400.10.20.30.40.50.6Normalized and Scaled Largest EigenvalueProbabilityβ=2−6 −4 −2 0 2 400.10.20.30.40.50.6Normalized and Scaled Largest EigenvalueProbabilityβ=4Figure 1: Probability distribution of scaled largest eigenvalue (105 repeti-tions, n = 109) with the boundary condition q(s)∼ Ai(s), as s → ∞. (6) The probability distribution f2(s), corresponding to β = 2, is then given by d f2(s) = ds F2(s), (7) where F2(s) = exp � − � ∞(x − s)q(x)2 dx � . (8) s The distributions for β = 1 and β = 4 are the derivatives of R∞sF1(s)2 = F2(s)e− q(x) dx (9) and �2R R� �21∞se−1 2∞sq(x) dxq(x) dx +2s e F2(s) (10)F4 = 2 . 22 3 To solve this numerically using MATLAB, first rewrite (5) as a first-order system: dq q�= ds q� sq + 2q3 . (11) 4� � � � � � � This can be solved as an initial-value problem starting at s = s0 = suf-ficiently large positive number, and integrating