An eigenvalue algorithm¶
$$ \newcommand{\eg}{{\it e.g.}} \newcommand{\ie}{{\it i.e.}} \newcommand{\argmin}{\operatornamewithlimits{argmin}} \newcommand{\mc}{\mathcal} \newcommand{\mb}{\mathbb} \newcommand{\mf}{\mathbf} \newcommand{\minimize}{{\text{minimize}}} \newcommand{\diag}{{\text{diag}}} \newcommand{\cond}{{\text{cond}}} \newcommand{\rank}{{\text{rank }}} \newcommand{\range}{{\mathcal{R}}} \newcommand{\null}{{\mathcal{N}}} \newcommand{\tr}{{\text{trace}}} \newcommand{\dom}{{\text{dom}}} \newcommand{\dist}{{\text{dist}}} \newcommand{\R}{\mathbf{R}} \newcommand{\SM}{\mathbf{S}} \newcommand{\ball}{\mathcal{B}} \newcommand{\bmat}[1]{\begin{bmatrix}#1\end{bmatrix}} \newcommand{\loss}{\ell} \newcommand{\eloss}{\mc{L}} \newcommand{\abs}[1]{| #1 |} \newcommand{\norm}[1]{\| #1 \|} \newcommand{\tp}{T} $$We explore an eigenvalue algorithm: given a diagonalizable matrix ${\displaystyle A}$, the algorithm will produce a number ${\displaystyle \lambda }$, which is the greatest (in absolute value) eigenvalue of ${\displaystyle A}$, and a nonzero vector ${\displaystyle v}$, which is a corresponding eigenvector of ${\displaystyle \lambda }$, that is, ${\displaystyle Av=\lambda v}$.
The algorithm starts with a vector ${\displaystyle q_{0}}$, which may be an approximation to the dominant eigenvector or a random vector. The method is described by the recurrence relation
$${\displaystyle q_{k+1}={\frac {A q_{k}}{\|A q_{k}\|}}}$$So, at every iteration, the vector ${\displaystyle q_{k}}$ is multiplied by the matrix ${\displaystyle A}$ and normalized.
If we assume ${\displaystyle A}$ has an eigenvalue that is strictly greater in magnitude than its other eigenvalues and the starting vector ${\displaystyle q_{0}}$ has a nonzero component in the direction of an eigenvector associated with the dominant eigenvalue, then a subsequence ${\displaystyle \left(q_{k}\right)}$ converges to an eigenvector associated with the dominant eigenvalue.
Under the two assumptions listed above, the sequence ${\displaystyle \left(\rho _{k}\right)}$ defined by
$${\displaystyle \rho_{k}={\frac {q_{k}^{T}A q_{k}}{q_{k}^{T}q_{k}}}}$$converges to the dominant eigenvalue.
import numpy as np
np.random.seed(3001)
A = np.random.randn(100,100)
A = A.T@A
(Problem 1) Define a function ase3001_ev_algorithm() that receives a symmetric matrix $A$ and returns the largest eigenvalue of $A$. Repeat the algorithm until,
Report the largest eigenvalue of $A$, and the number of iterations required. Check how rapidly $e_k$ diminishes in a $\log_{10}$ scale plot.
# your code here
(Problem 2) Use ase3001_ev_algorithm() to find all eigenvalues of $A$. Report all eigenvalues of $A$ and check if your answer is correct.
# your code here