#
# See the exercise "RandomMatrixTheory.pdf" from RandomMatrixTheory.html
# in http://www.physics.cornell.edu/~myers/teaching/ComputationalMethods/ComputerExercises/
#
"""Random Matrix Theory exercise"""
import scipy, scipy.linalg, pylab
def GOE(N):
"""
Creates an NxN element of the Gaussian Orthogonal Ensemble,
by creating an NxN matrix of Gaussian random variables
using the random array function scipy.random.standard_normal([shape])
with [shape] = (N,N).
and adding it to its transpose (applying scipy.transpose to the matrix)
#
Typing GOE(4) a few times, check it's symmetric and random
"""
m = scipy.random.standard_normal((N,N))
m = m + scipy.transpose(m)
return m
def GOE_Ensemble(num, N):
"""
Creates a list "ensemble" of length num of NxN GOE matrices.
Starts with ensemble equalling an empty list []
Then, for n in range(num) uses the function ensemble.append
to add GOE(N) to the list.
#
Check GOE_Ensemble(3,2) gives 3 2x2 symmetric matrices
"""
# ensemble = [GOE(N) for n in range(num)]
ensemble = []
for n in range(num):
ensemble.append(GOE(N))
return ensemble
def CenterEigenvalueDifferences(ensemble):
"""
For each symmetric matrix in the ensemble, calculates the difference
between the two center eigenvalues of the ensemble, and returns the
differences as an array, after chopping off any tiny imaginary parts.
#
Given an ensemble of symmetric matrices,
finds their size N
len(m) gives the number of rows in m
starts with an empty list of differences
for each m in ensemble
finds the eigenvalues
scipy.linalg.eigvals(m) gives the eigenvalues of m
sorts them
scipy.sort(e) sorts a list from smallest to largest
appends eigenvalue[N/2] - eigenvalue[N/2-1] to the list
(Note: python lists go from 0 ... N-1, so for N=2 this
gives the difference between the only two eigenvalues)
Converts the list of differences diff into an array.
Lists are python objects: they can be appended to but not multiplied.
Arrays are scipy objects that look the same as lists, but
can't be appended and can be multiplied.
Use d = scipy.array(d) to convert from a list to an array.
Takes the real part of the difference array (rounding errors)
For a scipy d, d.real is the real part
returns the (real part of the) array of differences
#
Check
ensemble = GOE_Ensemble(3,2)
CenterEigenvalueDifferences(ensemble)
gives three positive numbers, that look like eigenvalues of the 2x2 matrices
"""
# Size of matrix
N = len(ensemble[0])
diffs = []
for mat in ensemble:
eigenvalues = scipy.sort(scipy.linalg.eigvals(mat))
diffs.append(eigenvalues[N/2]-eigenvalues[N/2-1])
diffs = scipy.array(diffs)
diffs = diffs.real
return diffs
def CenterDiffHistogram(ensemble, bins=30, showPlot=True):
"""
Calculates the center eigenvalue difference of an ensemble of
NxN matrices, using GOE_ensemble and CenterEigenvalueDifferences.
Finds the averages diffAve of the differences.
scipy.sum(d) gives the total; len(d) gives the num in the ensemble
Calculates an array of normalized differences, dividing diff by diffAve
Plots a histogram of the normalized differences
pylab.hist(d, bins=bins, normed=1) prepares a plot of a histogram with
total area one (normalized to a probability distribution).
The notation def f(x, y=3) in a function definition gives a default
value of 3 for the optional input variable y. The same notation
in calling the function f(x, y=2) will overload the default value
of y with 2. Thus hist(d, bins=bins) overloads the default
value of bins in the histogram with the input variable "bins".
pylab.show() should be run (to display the graph) if showPlot==True.
Warning: pylab currently freezes the display, so you'll need to
close the graph before typing more Python commands.
#
Check that
ensemble = GOE_Ensemble(3,2)
CenterDiffHistogram(ensemble)
gives spikes at the three eigenvalue differences you found with
CenterEigenvalueDifferences(ensemble)
Then run with ensembles of thousands of matrices, with size N from 2-10
"""
diffs = CenterEigenvalueDifferences(ensemble)
meanDiff = scipy.sum(diffs)/len(diffs)
diffsNormalized = diffs/meanDiff
pylab.hist(diffsNormalized, bins=bins, normed=1)
if showPlot==True:
pylab.show()
def Wigner(s):
"""
Returns the Wigner surmise for the probability distribution rho(s)
for the eigenvalue differences in the GOE ensemble,
rho(s) = (pi s/2) exp(-pi s^2/4)
Python provides a number of complications:
s^2 is s**2
In Python (and in some other languages), always divide by floats
unless you want the value truncated (s/2.0 not s/2)
pi and exp are not Python functions, but scipy functions
The Wigner surmise is correct only for 2x2 GOE matrices, but it's
easily derived and is quite close to the true (much messier) answer
for larger matrices.
#
Test that Wigner(1.0) is reasonable, that Wigner(0.0) is zero,
and that Wigner(10.0) is small.
"""
return (scipy.pi*s/2.0) * scipy.exp(-scipy.pi*s**2/4.)
def CompareEnsembleWigner(ensemble,bins=30):
"""
Plots the center eigenvalue difference histogram using CenterDiffHistogram
with showPlot=False.
Then creates an array of s values using scipy.arange from 0.0 to 4.0
with a step of 0.01.
If the function Wigner is defined using only multiplication, powers,
and simple scipy functions like exponentials, Wigner(sArray) should
return an array of the Wigner distribution evaluated at the values of s.
pylab.plot(sArray,theory) should provide the theory curve;
pylab.show() then displays both histogram and Wigner surmise.
#
Test the agreement with thousands of matrices and N=2.
Test that the agreement remains good with N from 4-10, and at least 30 bins.
"""
CenterDiffHistogram(ensemble, bins=bins, showPlot=False)
sArray = scipy.arange(0.0,4.0,0.01)
theory = Wigner(sArray)
pylab.plot(sArray, theory)
pylab.show()
def PM1(N):
"""
Creates a symmetric NxN matrix with random entries of +-1.
I generated an asymmetric array scipy.random.randint, and then set the
bottom left equal to the top right in a double loop.
#
scipy.random.randint(min, max+1, (N,N)) will generate random integers
in the range (min, min+1, ..., max) inclusive.
(Notice that the range starts at min but ends one BEFORE max:
this makes sense in python, where length-N arrays start at zero
and end at N-1.)
You'll want to generate a random matrix of zeros and ones,
and then multiply by two and subtract one.
#
Symmetrizing needs i in range(N) and j in range(i)
[range, like arrays and randint, stops one before the end]
#
Test first, by generating some matrices and making sure they're symmetric
and +-1.
"""
m = scipy.random.randint(0,2,(N,N))*2.0-1.0
for i in range(N):
for j in range(i):
m[j,i] = m[i,j]
return m
def PM1_Ensemble(num, N):
"""
Generates a +-1 ensemble, as for GOE_Ensemble above.
Use with CompareEnsembleWigner to test that N=2 looks different
from the Wigner ensemble, but larger N looks close to the GOE ensemble.
#
This is a powerful truth: the eigenvalue differences of very general classes
of symmetric NxN matrices all have the same probability distribution
as N goes to infinity. This is what we call universality.
"""
# ensemble = [PM1(N) for n in range(num)]
ensemble = []
for n in range(num):
ensemble.append(PM1(N))
return ensemble
def yesno():
response = input(' Continue? (y/n) ')
if len(response)==0: # [CR] returns true
return True
elif response[0] == 'n' or response[0] == 'N':
return False
else: # Default
return True
def demo():
"""Demonstrates solution for exercise: example of usage"""
print("Random Matrix Theory Demo")
print(" GOE 2x2 vs. Wigner")
CompareEnsembleWigner(GOE_Ensemble(1000,2))
if not yesno(): return
print(" +-1 2x2 vs. Wigner")
CompareEnsembleWigner(PM1_Ensemble(1000,2))
if not yesno(): return
print(" +-1 10x10 vs. Wigner")
CompareEnsembleWigner(PM1_Ensemble(1000,10))
if __name__=="__main__":
demo()
# Copyright (C) Cornell University
# All rights reserved.
# Apache License, Version 2.0