Repeated matrix multiplication makes the column vectors converge to eigen vectors of the matrix

convergence_to_eigenvector.py

# ! pip install celluloid
import numpy as np
import matplotlib.pyplot as plt
from celluloid import Camera

def plot_mat(A, evals=None, evecs=None, fig=None, ax=None):
  if ax is None:
    fig, ax = plt.subplots(1, 2, figsize=(8, 4))
  A = A / np.linalg.norm(A, axis=0, keepdims=True)
  ax[0].imshow(A)
  ax[1].quiver([0, 0], [0, 0], A[0, :], A[1, :], angles='xy', scale_units='xy', color=["r", "b"], scale=1, alpha=0.5)
  axis_max = np.abs([A.min(), A.max()]).max()
  if evecs is not None:
    axis_max = np.abs([axis_max, evecs.min(), evecs.max()]).max()
    ax[1].quiver([0, 0], [0, 0], evecs[0, :], evecs[1, :], angles='xy', scale_units='xy', color="0.5", scale=1, alpha=0.5)
    ax[1].quiver([0, 0], [0, 0], -evecs[0, :], -evecs[1, :], angles='xy', scale_units='xy', color="0.5", scale=1, alpha=0.5)

  axis_max = min(axis_max*1.1, axis_max + 1)
  
  ax[1].set_xlim([-axis_max, axis_max])
  ax[1].set_ylim([-axis_max, axis_max])
  return fig, ax

if __name__ == "__main__":
  # rs = np.random.RandomState(42) # Good
  rs = np.random.RandomState(1201)
  A = rs.randn(2, 2)
  A = A @ A.T
  # A = A / np.linalg.norm(A, axis=0, keepdims=True)
  evals, evecs = np.linalg.eig(A)
  print(evals, evecs)
  plot_mat(A, evals=evals, evecs=evecs)

  fig, ax = plt.subplots(1, 2, figsize=(16, 8))
  _, _ = plot_mat(A, evals=evals, evecs=evecs, fig=fig, ax=ax)
  camera = Camera(fig)

  B = np.eye(A.shape[0])
  for i in range(20):
    B_new = B @ A
    # B_new = B_new / np.linalg.norm(B_new, axis=0, keepdims=True) # Keeps norm 1
    evals, evecs = np.linalg.eig(B_new)

    _, _ = plot_mat(B_new, evals=evals, evecs=evecs, fig=fig, ax=ax)  
    camera.snap()
    dist = np.linalg.norm(B - B_new)
    B = B_new
    # ax[1].set_title(f"i={i}, dist={dist:.3f}")
    ax[0].set_title(r"$B = B \times A$")
    ax[1].set_title(f"i={i}. Grey are eigen vectors.")
    if dist < 0.001:
      break

  animation = camera.animate()
  animation.save('animation.mp4', dpi=300, fps=2)