%% Eigenvalues of random operators
% Yuji Nakatsukasa, April 2017

%%
% (Chebfun example ode-eig/Randfuneig.m)

%% 1. Eigenvalues of random matrices
% As is well known [2,3], eigenvalues of large random matrices 
% are distributed uniformly in a disk.
% This is called the circular law and it holds very generally (see [4]); 
% the matrix elements are allowed to come from a variety of random
% variables, as long as they are independent (a property called universality). Here we verify
% the circular law with a standard Gaussian random matrix. 
% Throughout this example, we scale the entries 
% in the standard fashion so that the asymptotic spectral radius converges to 1. 

n = 1000;
A = randn(n)/sqrt(n);
plot(eig(A),'k.'),hold on
plot(chebfun(@(x)exp(1i*x),[0 2*pi]))
axis equal off

%% 2. Eigenvalues of random low-rank matrices
% Low-rank matrices are everywhere these days. A natural question arises:
% what are the eigenvalues of random low-rank matrices? 
% Specifically, consider matrices given by $AB^T$ where $A,B$ are
% tall skinny $m\times n$ random matrices. Of course there are $m-n$ eigenvalues at 0.
% Where are the other eigenvalues? An important fact here 
% (which becomes crucial below)
% is that the nonzero eigenvalues are the same as those of
% $B^TA$, which is $n\times n$, hence smaller. 

n = 1000;
m = n*10;   % aspect ratio 10
A = randn(m,n)/power(m*n,1/4); % scaled random rectangular matrix
B = randn(m,n)/power(m*n,1/4);

hold off
plot(eig(B'*A),'k.'), hold on
plot(chebfun(@(x)exp(1i*x),[0 2*pi]))
axis equal off

%%
% The eigenvalues again appear to be distributed uniformly on the unit disk! 
% A difficulty in explaining this is that here the elements aren't exactly independent. 
% Most results in the literature on random matrices concern i.i.d.
% elements. Nonetheless, recent studies [1,5] explain the phenomenon:
% the eigenvalues of $AB^T$ are indeed distributed in the unit disk, as
% $m,n\rightarrow \infty$,
% but generally not uniformly. 
% More precisely, denoting by $k=m/n$ the aspect ratio, 
% the eigenvalue density at radius $r(0<r<1)$ converges as $n\rightarrow \infty$
% almost surely to $g(r,k) = k/\sqrt{(1-k)^2+4kr^2}$. Note that as $k\rightarrow \infty$, 
% $g(r,k)$ converges to the
% constant $1$, independent of $r$. 
% This indicates that as we take the limit where the 
% aspect ratio tends to infinity (while $m,n\rightarrow \infty$), the
% eigenvalue distribution converges to the uniform distribution over the unit disk. 

%%
% Here is a qualitative explanation: As $m/n\rightarrow \infty$, the columns of 
% the matrices $A,B$ become close to orthonormal, so there is an
% orthogonal $m\times m$ matrix $Q$ such that
% $QA\approx \big[I_n, 0\big]^T$ ($Q$ would also be random),
% and then the eigenvalues $AB^T$, which are the same as those of $(QB)^T(QA)$,
% are approximately those of the top $n\times n$ block of $QB$, which has 
% elements that are i.i.d. 

%% 3. Eigenvalues of a product of two random matrices
% The function $g(r,k)$ above suggests that for a product of two random square matrices, 
% the eigenvalues should be markedly clustered near the origin. Let's
% illustrate this:

n = 1000;
m = n; % we just changed this from n*10 in the low-rank case
A = randn(m,n)/power(m*n,1/4); 
B = randn(m,n)/power(m*n,1/4);

hold off
plot(eig(B'*A),'k.'), hold on
plot(chebfun(@(x)exp(1i*x),[0 2*pi]))
axis equal off

%%
% One can play with the aspect ratio and see that the larger $m/n$ is, the
% more uniform the eigenvalue distribution becomes. 


%% 4. Eigenvalues of a Fredholm operator defined by a random, finite-rank bivariate function
% Now we explore the continuous analogue. 
% We construct a random bivariate function, generated by a new Chebfun
% command randnfun2: this generates a random bivariate function 
% obtained by forming a function of random Fourier coefficients (with support in a disk), 
% then restricting the domain to a smaller square. 
% Let's generate such a function and look at the associated eigenvalues: here,
% we examine the eigenvalues of the Fredholm integral operator whose kernel
% is the bivariate function (this is a compact operator).
% Here, the eig(chebfun2) code computes the nonzero eigenvalues, using again (only more crucially)
% the identity $eig(AB) = eig(BA)$, as there are infinitely many eigenvalues at 0). 

dt = 0.01;                     % max wavenumber 2pi/dt
dom = [-1 1 -1 1];
f = randnfun2(dom,dt,'norm');  % normalized randum function 
ei = eig(f);
disp(['Number of nonzero eigenvalues: ',num2str(length(ei))])

hold off, plot(ei,'k.'), hold on
plot(chebfun(@(x)exp(1i*x),[0 2*pi]))
axis equal off

%%
% Again we see eigenvalues in a disk; are they uniformly distributed? 
% By analogy with the finite-dimensional case, perhaps we expect the eigenvalues to be uniform, as we can regard $f(x,y)$ as $A(y)B(x)^T$ where $A,B$ are quasimatrices, i.e.,
% matrices of size $\infty \times n$, that is, the aspect ratio is
% $m/n=\infty$. 
% However, we think the distribution here is not quite uniform. 
% 
%%
% Experiments suggest that the uniformness depends on the support of the
% Fourier coefficients used for randnfun2. If we change it to bounded maximum
% degree, we obtain a more convincingly uniform distribution, aside from some clustering along the real axis: 
% (below we extract the essence of randnfun and change the support of the coefficients from a disk (2-norm bound) 
% to a square (max-norm bound); a 1-norm bound results in more clustered
% eigenvalues). It is an open problem to explain this observation! 

n = round(2/dt);
c = randn(2*n+1) + 1i*randn(2*n+1);     % random coefficients on a square
[x,y] = meshgrid(-n:n,-n:n);
if ( n>0 )
 %c = c.*((x/n).^2 + (y/n).^2 <= 1);    % 2-norm bound; default in randnfun 
 c = c.*(max(abs(x/n),abs(y/n)) <= 1);  % max-norm bound 
 %c = c.*( (abs(x/n)+abs(y/n)) <= 1);   % 1-norm bound: very clustered eigvals
end
c = c/sqrt(nnz(c));                     % ensure var = 1 at each point
n = round(1.2*(dom(4)-dom(3))/dt+2);
dom2 = [ dom(1) dom(1)+n*dt dom(3) ...  % bigger domain
    dom(3)+n*dt ]; 
f = chebfun2(c, dom2, 'coeffs', 'trig');% periodic chebfun2
f = real(f);
f = f{-1,1,-1,1}/sqrt(dt);              % restrict and normalize

ei = eig(f);
disp(['Number of nonzero eigenvalues: ',num2str(length(ei))])
hold off, plot(ei,'k.'), hold on
plot(chebfun(@(x)exp(1i*x),[0 2*pi]))
axis equal off

%% 5. References
% 
% 1. Burda, Z. and Jarosz, A. and Livan, G. and Nowak, M. A. and Swiech, A.,
% Eigenvalues and singular values of products of rectangular Gaussian
% random matrices,
% _Physical Review E._, 82(6),061-114,2010.
% 
% 2. J. Ginibre, Statistical ensembles of complex, quaternion, and real 
% matrices,
% _J. Math. Phys._,6, 440, 1965.
% 
% 3. M. L. Mehta,
% _Random Matrices_, Academic Press, 2004. 
%
% 4. T. Tao, V. Vu, Random matrices: Universality of ESDs and the circular law, 
% _The Annals of Probability_, 38(5), 2023--2065, 2010,
% 
% 5. X. Zeng, Limiting empirical distribution for eigenvalues of products of random rectangular matrices, 
% _Statistics & Probability Letters_, 126, 33--40, 2017.  
% 
