See lecture notes and documentation on Brightspace for Python and Jupyter basics. If you are stuck, try to google or get in touch via Discord.
+Solutions must be submitted via the Jupyter Hub.
+Make sure you fill in any place that says YOUR CODE HERE or "YOUR ANSWER HERE".
In the following you will implement your own eigenvalue / eigenvector calculation routines based on the inverse power method and the iterated QR decomposition.
+ +import numpy as np
+import math as m
+We start by implementing the inverse power method, which calculates the eigenvector corresponding to an eigenvalue which is closest to a given parameter $\sigma$. In detail, you should implement a Python function $\text{inversePower(A, sigma, eps)}$, which takes as input the $n \times n$ square matrix $A$, the parameter $\sigma$, as well as some accuracy $\varepsilon$. It should return the eigenvector $\mathbf{v}$ (corresponding to the eigenvalue which is closest to $\sigma$) and the number of needed iteration steps.
+To do so, implement the following algorithm. Start by setting up the needed input:
+\begin{align} + B &= \left( A - \sigma \mathbf{1} \right)^{-1} \\ + \mathbf{b}^{(0)} &= (1,1,1,...) +\end{align}where $\mathbf{b}^{(0)}$ is a vector with $n$ entries. Then repeat the following and increase $k$ each iteration until the error $e$ is smaller than $\varepsilon$:
+\begin{align} + \mathbf{b}^{(k)} &= B \cdot \mathbf{b}^{(k-1)} \\ + \mathbf{b}^{(k)} &= \frac{\mathbf{b}^{(k)}}{|\mathbf{b}^{(k)}|} \\ + e &= \sqrt{ \sum_{i=0}^n \left(|b_i^{(k-1)}| - |b_i^{(k)}|\right)^2 } +\end{align}Return the last vector $\mathbf{b}^{(k)}$ and the number of needed iterations $k$.
+Verify your implementation by calculating all the eigenvectors of the matrix below and comparing them to the ones calculated by $\text{numpy.linalg.eig()}$. Then implement a unit test for your function.
+\begin{align*} + A = \begin{pmatrix} + 3 & 2 & 1\\ + 2 & 3 & 2\\ + 1 & 2 & 3 + \end{pmatrix}. +\end{align*} +def inversePower(A, sigma, eps):
+ """
+ Estimates the eigenvectors of matrix A by the inverse power method.
+
+ Args:
+ A: an n x n matrix
+ sigma: an initial guess for an eigenvector
+ eps: the desired accuracy
+
+ Returns:
+ A tuple (b, k) is returned, containing:
+ b: the eigenvector b corresponding to the eigenvalue
+ closests to sigma after k iterations
+ k: the number of iterations done
+
+ See also:
+ https://www.sciencedirect.com/topics/mathematics/inverse-power-method
+ """
+
+ # Does https://johnfoster.pge.utexas.edu/numerical-methods-book/LinearAlgebra_EigenProblem1.html help?
+
+ # A should be n x n.
+ n = len(A)
+ assert len(A.shape) == 2 and A.shape[0] == A.shape[1]
+
+ # Setup some initial values.
+ #B = np.linalg.inv(A - sigma*np.ones(n))
+ B = np.linalg.inv(A - sigma*np.eye(n))
+ #B = 1/(A - sigma*np.ones(n))
+ #B = 1/(A - sigma*np.eye(n))
+ b = np.ones(n)
+ k = 0
+ e = 0
+
+ while e > eps or k == 0:
+ b_prev = b.copy()
+ k += 1
+
+ b = B @ b
+ b /= np.sqrt(b @ b) # although b = np.linalg.norm(b) could be used
+ e = np.sqrt(np.sum( (np.abs(b_prev) - np.abs(b))**2) )
+
+ return b, k
+# Use this cell for your own testing ...
+
+A = np.array([[3, 2, 1], [2, 3, 2], [1, 2, 3]])
+
+sigma_list = [.03, 6., 1.99999989999999]
+#sigma = 0.03
+
+for sigma in sigma_list:
+ print("For sigma =", sigma)
+ b, k = inversePower(A, sigma, 1e-6)
+ print("b =", b)
+ lam = np.dot(A, b) / b
+ print("lam =", lam)
+ print()
+ print()
+
+print(np.linalg.eig(A)[0])
+def test_inversePower():
+ A = np.array([[3, 2, 1], [2, 3, 2], [1, 2, 3]])
+
+test_inversePower()
+Next you will need to implement the tri-diagonalization scheme following Householder. To this end implement a Python function $\text{tridiagonalize(A)}$ which takes a symmetric matrix $A$ as input and returns a tridiagonal matrix $T$ of the same dimension. Therefore, your algorithm should execute the following steps:
+First use an assertion to make sure $A$ is symmetric. Then let $k$ run from $0$ to $n-1$ and repeat the following:
+\begin{align} + q &= \sqrt{ \sum_{j=k+1}^n \left(A_{j,k}\right)^2 } \\ + \alpha &= -\operatorname{sgn}(A_{k+1,k}) \cdot q \\ + r &= \sqrt{ \frac{ \alpha^2 - A_{k+1,k} \cdot \alpha }{2} } \\ + \mathbf{v} &= \mathbf{0} \qquad \text{... vector of dimension } n \\ + v_{k+1} &= \frac{A_{k+1,k} - \alpha}{2r} \\ + v_{k+j} &= \frac{A_{k+j,k}}{2r} \quad \text{for } j=2,3,\dots,n \\ + P &= \mathbf{1} - 2 \mathbf{v}\mathbf{v}^T \\ + A &= P \cdot A \cdot P +\end{align}At the end return $A$ as $T$.
+Apply your routine to the matrix $A$ defined in task 1 as well as to a few random (but symmetric) matrices of different dimension $n$.
+Hint: Use $\text{np.outer()}$ to calculate the matrix $\mathbf{v}\mathbf{v}^T$ needed in the definition of the Housholder transformation matrix $P$.
+ +def tridiagonalize(A):
+ """
+ Tridiagonalizes a given matrix following the Householder method.
+ Args:
+ A: NxN matrix
+ Returns:
+ The matrix A, but now tridiagonalized
+ """
+ assert A.shape[0] == A.shape[1] and len(A.shape) == 2
+ n = len(A)
+ for k in range(n-1):
+ q = m.sqrt(np.sum(A[k+1:n,k]**2))
+ alpha = -1*np.sign(A[k+1,k])*q
+ r = m.sqrt((alpha**2 - A[k+1,k]*alpha)/2)
+ v = np.zeros(n)
+ v[k+1] = (A[k+1,k]-alpha)/(2*r)
+ print(np.shape(v),np.shape(A),np.shape(r),np.shape(A[k+2:n]/(2*r)))
+ v[k+2:n] = A[k+2:n,k]/(2*r)
+ P = np.identity(n)-(2*np.outer(v,v))
+ A = np.dot(np.dot(P,A),P)
+ return A
+B = np.array([[3,2,1],[2,3,2],[1,2,3]])
+
+tridiagonalize(B)
+Implement the $QR$ decomposition based diagonalization routine for tri-diagonal matrices $T$ in Python as a function $\text{QREig(T, eps)}$. It should take a tri-diagonal matrix $T$ and some accuracy $\varepsilon$ as inputs and should return all eigenvalues in the form of a vector $\mathbf{d}$. By making use of the $QR$ decomposition as implemented in Numpy's $\text{numpy.linalg.qr()}$ the algorithm is very simple and reads:
+Repeat the following until the error $e$ is smaller than $\varepsilon$:
+\begin{align} + Q \cdot R &= T^{(k)} \qquad \text{... do this decomposition with the help of Numpy!} \\ + T^{(k+1)} &= R \cdot Q \\ + e &= |\mathbf{d_1}| +\end{align}where $T^{(0)} = T$ and $\mathbf{d_1}$ is the first sub-diagonal of $T^{(k+1)}$ at each iteration step $k$. At the end return the main-diagonal of $T^{(k+1)}$ as $\mathbf{d}$.
+Implement a unit test for your function based on the matrix $A$ defined in task 1. You will need to tri-diagonalize it first.
+ +def QREig(T, eps):
+ """
+ Follows the method of the QR decomposition based diagonalization routine for tridiagonalized matrices. The matrix T is
+ diagonalized, resulting in all diagonal elements being an eigenvalue.
+ Args:
+ T: a already tridiagonalized matrix
+ eps: the desired accuracy
+ Returns:
+ A one dimensional array with the eigenvalues of the matrix T
+ """
+ e = eps + 1
+ while e > eps:
+ Q,R = np.linalg.qr(T)
+ T = np.matmul(R,Q)
+ e = np.sum(np.absolute(np.diag(T,k=1)))
+ #print(np.diag(T))
+ return np.diag(T)
+A_tridiag = tridiagonalize(A)
+print(A_tridiag)
+def test_QREig():
+ """
+ Tests the QREig function for the matrix A_tridiag (defined in task 2) and eps = 0.001
+ """
+ eps = 0.001
+ x = QREig(A_tridiag,eps)
+ print(x)
+
+test_QREig()
+With the help of $\text{QREig(T, eps)}$ you can now calculate all eigenvalues. Then you can calculate all corresponding eigenvectors with the help of $\text{inversePower(A, sigma, eps)}$, by setting $\sigma$ to approximately the eigenvalues you found (you should add some small random noise to $\sigma$ in order to avoid singularity issues in the inversion needed for the inverse power method).
+Apply this combination of functions to calculate all eigenvalues and eigenvectors of the matrix $A$ defined in task 1.
+ +"""
+Eigenvalues and -vectors calculated numerically for matrix A using functions from previous tasks.
+"""
+eps = 0.001
+T = QREig(A_tridiag,eps)
+for i in range(len(T)):
+ sigma = T[i] + np.random.normal(0,0.01,1)
+ x = inversePower(A_tridiag,sigma,eps)
+ print(x)
+Test your eigenvalue / eigenvector algorithm for larger random (but symmetric) matrices.
+ +"""
+Completely 'random' and totally not self-made matrices C and D are being tridiagonalized and its eigenvalues and -vectors are
+calculated.
+"""
+C = np.array([[5,4,3,2,1]
+ ,[3,4,5,1,2]
+ ,[5,1,4,2,3]
+ ,[3,2,1,5,1]
+ ,[5,4,3,2,1]])
+#print(C)
+C_tridiag = tridiagonalize(C)
+T = QREig(C_tridiag,eps)
+for i in range(len(T)):
+ sigma = T[i] + np.random.normal(0,0.01,1)
+ x = inversePower(C_tridiag,sigma,eps)
+ print(x)
+
+D = np.array([[6,4,3,2,1]
+ ,[3,3,5,1,2]
+ ,[5,1,4,2,2]
+ ,[3,1,1,5,6]
+ ,[1,4,7,2,3]])
+#print(C)
+C_tridiag = tridiagonalize(C)
+T = QREig(C_tridiag,eps)
+for i in range(len(T)):
+ sigma = T[i] + np.random.normal(0,0.01,1)
+ x = inversePower(C_tridiag,sigma,eps)
+ print(x)
+