In [1]:
%matplotlib inline

Rayleigh quotient

  • $M$: Hermitian matrix

  • $x$: nonzero vector

  • Rayleigh quotient: $R(M,x)={x^{*}Mx \over x^{*}x}$

If $x$ is an eigenvector of $A$ then $R(M,x)$ is its eigenvalue. Otherwise, $R(M,x)$ is the closest scalar to an eigenvalue (in least squares):

Exercise 1:

Which is the value $\alpha$ that minimizes $||Ax - \alpha x||_2^2$?

hint: Compute the derivative of $||Ax - \alpha x||_2^2$ w.r.t. $\alpha$

Exercise 2:

  • Set numpy's random generator to 0 (use the np.random.seed function)
  • Generate $A$, a $2 \times 2$ random matrix (use np.random.randn)
  • Compute $M = A^* A$. Is this matrix Hermitian?
  • Compute its eigenvalues and eigenvectors (use np.linalg.eig)
  • Choose one of the eigenvectors and compute its Rayleigh quotient
  • Write a function that compute the Rayleigh quotient of a vector, w.r.t. a matrix, both sent as parameters
  • Generate 1000 random vectors $\sim \mathcal{N}(\mu = 0,\sigma = 10)$
  • For each vector compute its Rayleigh quotient, w.r.t. $M$
  • Plot the histogram of the resulting Rayleigh quotients (use matplotlib.pyplot.hist)
In [ ]:

Power iteration

A well known method to compute eigenvalues is the power iteration method (also known as the Von Mises iteration).

This method takes $A$ a diagonalizable matrix as an input, and it outputs its greatest eigenvalue $\lambda$ (in absolute value), and its corresponding eigenvector $v$ (i.e., $Av=\lambda v$).

Power iteration($A$):

  • initialization: create $b$ a random vector s.t. $||b||_2 = 1$
  • Repeat:
    • $w \leftarrow Ab$
    • $b \leftarrow w/||w||$
    • $\lambda = b^T A b$

This method does not converge if:

  • The initial vector $b$ is orthogonal with the eigenvector of $\lambda$
  • The eigenvalue is not unique

Exercise 3 How does this method work?

  • Write $b^{(k)}$ ($b$ at iteration $k$) as a function of $b^{(0)}$ ($b$ at iteration 0) and $A$
  • Let us assume that the initial vector is s.t.: $b^{0}=c_{1}v_{1}+c_{2}v_{2}+\cdots +c_{m}v_{m}$, with $v_i$ the $i$-th eigenvector.
  • Decompose $A^k b^{(0)}$ according to the previous equation.
  • Factorize by $\lambda_1$, i.e., the largest eigenvalue.
  • When $k$ is large to which value does the previous expression converge?

Exercise 4

  • Write the Power iteration method
  • In order to analyze the algorithm compute:
    • The eigenvalues and eigenvectors using np.linalg.eig
    • Compute the difference between the two largest eigenvalues
    • At each iteration compute the difference between the largest eigenvalues output by np.linalg.eig, and the current value estimated by the power iteration method.
  • For different random matrices $M$ generated as in the previous exercise:
    • Run 10 iterations of the power iteration method
    • Record and plot the evolution of the error
In [ ]:

Inverse iteration method

  • $A$: square matrix, $\forall i, \;\lambda_i$ are eigenvalues of $A$
  • Let $\mu \in \mathbb{R}$ s.t. $\mu \neq \lambda_i$

Matrix $(A − \mu I)^{−1}$ and $A$ have the same eigenvectors, and the eigenvalues of $(A − \mu I)^{−1}$ are $(\lambda - \mu)^{-1}$

Exercise 5: Prove the previous statement.

If $\mu$ is close to the eigenvalue $\lambda_j$ then $(\lambda_j - \mu \mathbf{I})^{-1}$ is a larger eigenvalue than $(\lambda_i - \mu \mathbf{I})^{-1}, \; \forall j \neq i$.

The power iteration method applied to $(\lambda_i - \mu \mathbf{I})^{-1}$ converges quickly.

Inverse Iteration Algorithm ($A$, $\mu$):

  • Initialization: Generate $v$ a random vector with norm 1
  • Repeat:
    • $w \leftarrow (A-\mu I)^{-1}v$
    • $v \leftarrow w/||w||$
    • $\lambda \leftarrow v^TAv$

Exercise 6:

  • Program the inverse iteration algorithm
  • Generate a random matrix $M$ like in the previous exercises
  • Test the algorithm setting $\mu = 10$ and iterating for 100 steps, and compare to the output of the np.linalg.eig function.
  • For 100 random values of $\mu ~\mathcal{N}(0,\sigma=10)$, compute an eigenvalue with the inverse iteration algorithm.
  • Print the eigenvalues found with 2 decimals precision (use the np.round function)
In [ ]:

Iteration of Rayleigh quotient

Instead of using a random value $\mu$ as initial guess, we can use the Rayleigh quotient instead.

Moreover updating $\mu$ at each iteration allows a faster convergence.

These ideas give birth to the Rayleigh quotient iteration method

Rayleigh quotient iteration($A$):

  • Initialization: Generate $v$ a random vector with norm 1
  • $\lambda \leftarrow v^T A v$ (Rayleigh quotient)
  • Repeat:
    • $w = (A-\lambda)^{-1}v$
    • $v = w/||w||$
    • $\lambda = v^T A v$

Exercise 7:

  • Program the Rayleigh quotient iteration algorithm
  • Generate a random matrix $M$ like in the previous exercises
  • Test the algorithm, and compare to the output of the np.linalg.eig function.
  • Compute 100 eigenvalues with the inverse iteration algorithm.
  • Print the eigenvalues found with 2 decimals precision (use the np.round function)
In [ ]:

Schur factorization

Goal: Factorize $A$ in order to reveal the eigenvalues.

$$A = Q T Q^*$$
  • $Q$: unitary matrix
  • $T$: Upper triangular matrix
  • $T$ and $A$ are similar matrices, so they have the same eigenvalues

Theorem: Every square matrix has a Schur decomposition

Method: This factorization is computed iteratively using the $QR$ factorization

Schur factorization(A):

  • Let $A_0 = A$
  • Repeat for k=0 to $NbIterations$:
    • Compute the factorization $A_k = Q_k R_k$
    • Define $A_{k+1} = R_k Q_k = Q_k^* A_k Q_k$

Exercise 8:

  • Program the Schur factorization algorithm
  • Generate a matrix $M$ like in the previous exercises/.
  • Compute the matrix $T$ and compare the diagonal of $T$ with the eigenvalues computed with the np.linalg.eig function
In [ ]: