From b6f0b0d7d94859e777936ff79a47fa5627092ac3 Mon Sep 17 00:00:00 2001 From: Kees van Kempen Date: Thu, 24 Feb 2022 18:03:36 +0100 Subject: [PATCH] Fetch week 4 --- Week 4/8 Eigenvalues and Eigenvectors.ipynb | 528 ++++++++++++++++++++ 1 file changed, 528 insertions(+) create mode 100644 Week 4/8 Eigenvalues and Eigenvectors.ipynb diff --git a/Week 4/8 Eigenvalues and Eigenvectors.ipynb b/Week 4/8 Eigenvalues and Eigenvectors.ipynb new file mode 100644 index 0000000..d810302 --- /dev/null +++ b/Week 4/8 Eigenvalues and Eigenvectors.ipynb @@ -0,0 +1,528 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "4ec40081b048ce2f34f3f4fedbb0be10", + "grade": false, + "grade_id": "cell-98f724ece1aacb67", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "# CDS: Numerical Methods Assignments\n", + "\n", + "- 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.\n", + "\n", + "- Solutions must be submitted via the Jupyter Hub.\n", + "\n", + "- Make sure you fill in any place that says `YOUR CODE HERE` or \"YOUR ANSWER HERE\".\n", + "\n", + "## Submission\n", + "\n", + "1. Name all team members in the the cell below\n", + "2. make sure everything runs as expected\n", + "3. **restart the kernel** (in the menubar, select Kernel$\\rightarrow$Restart)\n", + "4. **run all cells** (in the menubar, select Cell$\\rightarrow$Run All)\n", + "5. Check all outputs (Out[\\*]) for errors and **resolve them if necessary**\n", + "6. submit your solutions **in time (before the deadline)**" + ] + }, + { + "cell_type": "raw", + "metadata": {}, + "source": [ + "team_members = \"\"" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "26a1ebc121c0cc223b0079b3fcb1d606", + "grade": false, + "grade_id": "cell-79ad6e3d4ad2eb15", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "## Eigenvalues and Eigenvectors\n", + "\n", + "In the following you will implement your own eigenvalue / eigenvector calculation routines based on the inverse power method and the iterated QR decomposition." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "5727ddddb338b9fd153d42af584e4e92", + "grade": true, + "grade_id": "cell-c6c3a2556ae5174d", + "locked": false, + "points": 0, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "# Import packages here ...\n", + "\n", + "# YOUR CODE HERE\n", + "raise NotImplementedError()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "1bcdeef3b38d48f61ea06b7c63fed895", + "grade": false, + "grade_id": "cell-ed71eb4882ee221a", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "### Task 1\n", + "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.\n", + "\t\n", + "To do so, implement the following algorithm. Start by setting up the needed input:\n", + "\n", + "\\begin{align}\n", + " B &= \\left( A - \\sigma \\mathbf{1} \\right)^{-1} \\\\\n", + " \\mathbf{b}^{(0)} &= (1,1,1,...)\n", + "\\end{align}\n", + "\n", + "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$:\n", + "\n", + "\\begin{align}\n", + " \\mathbf{b}^{(k)} &= B \\cdot \\mathbf{b}^{(k-1)} \\\\\n", + " \\mathbf{b}^{(k)} &= \\frac{\\mathbf{b}^{(k)}}{|\\mathbf{b}^{(k)}|} \\\\\n", + " e &= \\sqrt{ \\sum_{i=0}^n \\left(|b_i^{(k-1)}| - |b_i^{(k)}|\\right)^2 }\n", + "\\end{align}\n", + "\n", + "Return the last vector $\\mathbf{b}^{(k)}$ and the number of needed iterations $k$. \n", + "\n", + "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.\n", + "\n", + "\\begin{align*}\n", + " A = \\begin{pmatrix}\n", + " 3 & 2 & 1\\\\ \n", + " 2 & 3 & 2\\\\\n", + " 1 & 2 & 3\n", + " \\end{pmatrix}.\n", + "\\end{align*}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "fbd7848114b10e527cfdadae70432b75", + "grade": true, + "grade_id": "cell-876c5ab872a98046", + "locked": false, + "points": 3, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "def inversePower(A, sigma, eps):\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "d12ec4311abc0ec7ad2f4e7c12b6a247", + "grade": true, + "grade_id": "cell-e12bd2509520938e", + "locked": false, + "points": 0, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "# Use this cell for your own testing ...\n", + "\n", + "# YOUR CODE HERE\n", + "raise NotImplementedError()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "4848b41fcb6770679285941320a2de3b", + "grade": true, + "grade_id": "cell-f458054aeff75cc9", + "locked": false, + "points": 3, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "def test_inversePower():\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()\n", + " \n", + "test_inversePower()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "4c2230bb7698096ad982513c84b1b013", + "grade": false, + "grade_id": "cell-d25e9b53cd92f821", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "### Task 2 \n", + "\n", + "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:\n", + "\t\n", + "First use an assertion to make sure $A$ is symmetric. Then let $k$ run from $0$ to $n-1$ and repeat the following:\n", + "\n", + "\\begin{align}\n", + " q &= \\sqrt{ \\sum_{j=k+1}^n \\left(A_{j,k}\\right)^2 } \\\\\n", + " \\alpha &= -\\operatorname{sgn}(A_{k+1,k}) \\cdot q \\\\\n", + " r &= \\sqrt{ \\frac{ \\alpha^2 - A_{k+1,k} \\cdot \\alpha }{2} } \\\\\n", + " \\mathbf{v} &= \\mathbf{0} \\qquad \\text{... vector of dimension } n \\\\\n", + " v_{k+1} &= \\frac{A_{k+1,k} - \\alpha}{2r} \\\\\n", + " v_{k+j} &= \\frac{A_{k+j,k}}{2r} \\quad \\text{for } j=2,3,\\dots,n \\\\\n", + " P &= \\mathbf{1} - 2 \\mathbf{v}\\mathbf{v}^T \\\\\n", + " A &= P \\cdot A \\cdot P\n", + "\\end{align}\n", + "\n", + "At the end return $A$ as $T$.\n", + "\n", + "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$.\n", + "\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$. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "5d057c8a5c00f59b856d55de8b2e5923", + "grade": true, + "grade_id": "cell-0cbf043c8cbc0601", + "locked": false, + "points": 3, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "def tridiagonalize(A):\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "3eab15fa13ffbb2cd9383870f305e206", + "grade": true, + "grade_id": "cell-460dfaeef80dd101", + "locked": false, + "points": 0, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "# Apply your routine here ...\n", + "\n", + "# YOUR CODE HERE\n", + "raise NotImplementedError()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "b847d4c4fa2311288742cf0badc5a053", + "grade": false, + "grade_id": "cell-e1a2cdfe8e8c0bd6", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "### Task 3\n", + "\n", + "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:\n", + "\n", + "Repeat the following until the error $e$ is smaller than $\\varepsilon$:\n", + "\n", + "\\begin{align}\n", + " Q \\cdot R &= T^{(k)} \\qquad \\text{... do this decomposition with the help of Numpy!} \\\\\n", + " T^{(k+1)} &= R \\cdot Q \\\\\n", + " e &= |\\mathbf{d_1}| \n", + "\\end{align}\n", + "\n", + "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}$. \n", + "\n", + "Implement a unit test for your function based on the matrix $A$ defined in task 1. You will need to tri-diagonalize it first." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "bb2814035b4f1e593517cb2a54411b38", + "grade": true, + "grade_id": "cell-17fe57c93a9ad465", + "locked": false, + "points": 3, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "def QREig(T, eps):\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "a2ce6f6cdaa5b7e83b5bdb2fe3b12370", + "grade": true, + "grade_id": "cell-19e1c89cfea6b4ca", + "locked": false, + "points": 0, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "# Use this cell for your own testing ...\n", + "\n", + "# YOUR CODE HERE\n", + "raise NotImplementedError()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "0bc7766605e0934f89b9df038f82dc51", + "grade": true, + "grade_id": "cell-e27978374106d1a4", + "locked": false, + "points": 3, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "def test_QREig():\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()\n", + " \n", + "test_QREig()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "f1725d3ecdd5358e4129b9006a4b8bc1", + "grade": false, + "grade_id": "cell-c530435b67593096", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "### Task 4\n", + "\n", + "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). \n", + "\n", + "Apply this combination of functions to calculate all eigenvalues and eigenvectors of the matrix $A$ defined in task 1." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "8792c616fbe2b5159caed5e873050243", + "grade": true, + "grade_id": "cell-bb8315bb7895761f", + "locked": false, + "points": 3, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "# YOUR CODE HERE\n", + "raise NotImplementedError()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "1ea6ea8324c01aee6a8e3225cbcd98a5", + "grade": false, + "grade_id": "cell-a1926d9a9975679b", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "### Task 5\n", + "\n", + "Test your eigenvalue / eigenvector algorithm for larger random (but symmetric) matrices." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "42126fc8bf851a1b021f8ecf7cb2bc6f", + "grade": true, + "grade_id": "cell-9be1ee8fbf7ec7a8", + "locked": false, + "points": 2, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "# YOUR CODE HERE\n", + "raise NotImplementedError()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.10" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}