import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation, rc
from IPython.display import HTML
import warnings
warnings.filterwarnings("ignore", category=np.ComplexWarning)
import sys
sys.path.append('../')
import RandomMatrix as RM
np.random.seed(140)
We will now create a matrix-valued stochastic process, specifically a GUE process. Define $A_t$ as the matrix state of the process at time $t$. If $0 \leq t_1 < t_2 < t_3$, then $A_{t_2} - A_{t_1} \stackrel{d}{=} \sqrt{t_2-t_1} B$ where $B$ is drawn from the GUE and $A_{t_2} - A_{t_1}$ is independent from $A_{t_3} - A_{t_2}$. We will choose our initial condition ($A_0$) to be the zero matrix.
Next we will choose $t_{\text{max}}$ such that we can run our process from $0 \leq t \leq t_{\text{max}}$
t_max = 100
Next we generate a 2x2 symmetric standard normal matrices that will represent the displacements of our process through time. We will measure increments after one unit of time so the increments will simply be from a generated GUE.
init_matrix = np.zeros((2,2))
matrix_displacements = [RM.Generate_GUE(2) for time in range(t_max-1)]
Similar to regular Brownian Motion, the process can be created by taking the cumulative sum of our displacement values.
process_data = [init_matrix]
process_data.extend(matrix_displacements)
A_t = np.cumsum(process_data, axis=0)
We can now calculate the eigenvalues of the process at every time.
A_eigs = np.linalg.eigvals(A_t)
A_eigs.sort()
We can visualize how to eigenvalues behave on the real axis (recall we are currently working with Hermitian matrices, so the eigenvalues must be real).
fig, ax = plt.subplots()
ax.set_xlim(-50,50)
ax.set_ylim(-1,1)
plt.xlabel("Re")
plt.ylabel("Im")
plt.title("Process of Eigenvalues")
plt.plot([-50, 50], [0, 0], linewidth=1, c="black")
eig1, = ax.plot(A_eigs[0][0],0, marker="o")
eig2, = ax.plot(A_eigs[0][1],0, marker="o")
def update(t):
eig1.set_data([A_eigs[t][0]],[0])
eig2.set_data([A_eigs[t][1]],[0])
return (eig1,eig2)
eig_animation = animation.FuncAnimation(fig, update, interval=50, blit=True, repeat=True,
frames=100)
plt.close()
HTML(eig_animation.to_html5_video())
We can visualize a line plot of the movement of these eigenvalues parameterized by time.
fig, ax = plt.subplots()
ax.set_xlim(0,99)
ax.set_ylim(-50,50)
line1, = ax.plot([], [], lw=2)
line2, = ax.plot([], [], lw=2)
plt.xlabel("t")
plt.ylabel("Eigenvalues")
plt.title("Process of Eigenvalues")
eig1 = [max(B) for B in A_eigs]
eig2 = [min(B) for B in A_eigs]
def eig_path(time):
line1.set_data(np.arange(time), eig1[:time])
line2.set_data(np.arange(time), eig2[:time])
return (line1, line2)
eig_animation = animation.FuncAnimation(fig, eig_path, interval=50, blit=True, repeat=True,
frames=100)
plt.close()
HTML(eig_animation.to_html5_video())
The motion of these eigenvalues is known as Dyson Brownian Motion. The process of the path satisfies the following Stochastic Differential Equation
$$d\lambda_i = dB_i + \sum_{1 \leq j \leq n: j \neq i} \frac{dt}{\lambda_i - \lambda_j} $$We can run the above process for matrices of any size. Below we plot the movement of eigenvalues of the process for 5x5 matrices.
n = 5
init_5 = np.zeros((5,5))
A_5_displacements = [RM.Generate_GUE(n) for time in range(t_max-1)]
A_5_process = [init_5]
A_5_process.extend(A_5_displacements)
A_5_path = np.cumsum(A_5_process, axis=0)
A_5_eigs = np.linalg.eigvals(A_5_path)
A_5_eigs.sort()
fig, ax = plt.subplots()
ax.set_xlim(0,99)
ax.set_ylim(-50,50)
plt.xlabel("t")
plt.ylabel("Eigenvalues")
plt.title("Process of Eigenvalues")
lines = [ax.plot([], [], lw=2)[0] for i in range(n)]
def eig_path(time):
for i in range(n):
current_line = lines[i]
current_line.set_data(np.arange(time), A_5_eigs.T[i][:time])
return lines
eig_n_animation = animation.FuncAnimation(fig, eig_path, interval=50, blit=True, repeat=True,
frames=100)
plt.close()
HTML(eig_n_animation.to_html5_video())