{ "cells": [ { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "# Problem 1 (20 pts)\n", "\n", "- (5 pts) Prove that $\\mathrm{vec}(AXB) = (B^\\top \\otimes A)\\, \\mathrm{vec}(X)$ if $\\mathrm{vec}(X)$ is a columnwise reshape of a matrix into a long vector. What does it change if the reshape is rowwise? **Note:** to make a columnwise reshape in Python one should use ```np.reshape(X, order='f')```, where the string ```'f'``` stands for the Fortran ordering.\n", "\n", "\n", "- (5 pts) Let matrices $A$ and $B$ have eigendecompositions $A = S_A\\Lambda_A S_A^{-1}$ and $B = S_B\\Lambda_B S^{-1}_B$. Find eigenvectors and eigenvalues of the matrix $A\\otimes I + I \\otimes B$.\n", "\n", "\n", "- (10 pts) Let $A = \\mathrm{diag}\\left(\\frac{1}{1000},\\frac{1}{999},\\dots \\frac{1}{2}, 1, 1000 \\right)$. Estimate analytically the number of iterations required to solve linear system with $A$ with the relative accuracy $10^{-4}$ using\n", " - Richardson iteration with the optimal choice of parameter\n", " - Chebyshev iteration\n", " - Conjugate gradient method" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Problem 2 (30 pts)\n", "\n", "## Spectral graph partitioning and inverse iteration\n", "\n", "\n", "Given connected graph $G$ and its corresponding graph Laplacian matrix $L = D - A$ with eigenvalues $0=\\lambda_1, \\lambda_2, ..., \\lambda_n$, where $D$ is its degree matrix and $A$ is its adjacency matrix, *Fiedler vector* is an eignevector correspondng to the second smallest eigenvalue $\\lambda_2$ of $L$. Fiedler vector can be used for graph partitioning: positive values correspond to one part of a graph and negative values to another.\n", "\n", "\n", "#### Inverse power method (15 pts)\n", "\n", "To find the Fiedler vector we will use the inverse iteration with adaptive shifts (Rayleigh quotient iteration). \n", "\n", "* Write down the orthoprojection matrix on the space orthogonal to the eigenvector of $L$, corresponding to the eigenvalue $0$ and prove (analytically) that it is indeed an orthoprojection.\n", " \n", "* Implement the spectral partitioning as the function ```partition```:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# INPUT:\n", "# A - adjacency matrix (scipy.sparse.csr_matrix)\n", "# num_iter_fix - number of iterations with fixed shift (int)\n", "# shift - (float number)\n", "# num_iter_adapt - number of iterations with adaptive shift (int) -- Rayleigh quotient iteration steps\n", "# x0 - initial guess (1D numpy.ndarray)\n", "# OUTPUT:\n", "# x - normalized Fiedler vector (1D numpy.ndarray)\n", "# eigs - eigenvalue estimations at each step (1D numpy.ndarray)\n", "# eps - relative tolerance (float)\n", "def partition(A, shift, num_iter_fix, num_iter_adapt, x0, eps): # 10 pts\n", " x = x0\n", " eigs = np.array([0])\n", " return x, eigs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Algorithm must halt before `num_iter_fix + num_iter_adapt` iterations if the following condition is satisfied $\\|\\lambda_k - \\lambda_{k-1}\\|_2 / \\|\\lambda_k\\|_2 \\leq \\epsilon$ at some step $k$.\n", "\n", "Do not forget to use the orthogonal projection from above in the iterative process to get the correct eigenvector.\n", "It is also a good idea to use ```shift=0``` before the adaptive stragy is used. This, however, is not possible since the matrix $L$ is singular, and sparse decompositions in ```scipy``` do not work in this case. Therefore, we first use a very small shift instead.\n", "\n", "* Generate a random `lollipop_graph` using `networkx` library and find its partition. [Draw](https://networkx.github.io/documentation/networkx-1.9/examples/drawing/labels_and_colors.html) this graph with vertices colored according to the partition.\n", "\n", "* Start the method with a random initial guess ```x0```, set ```num_iter_fix=0``` and comment why the method can converge to a wrong eigenvalue.\n", "\n", "#### Stanford Bunny (15 pts)\n", "\n", "Let us now find a partition of a large graph, obtained from the 3D model from [The Stanford 3D Scanning Repository](http://graphics.stanford.edu/data/3Dscanrep/):\n", "\n", "* Install the [trimesh](https://pypi.python.org/pypi/trimesh) library. In order to visualize 3D model smoothly you will also need to install `pyglet` library (see trimesh instructions). Alternatively, you can draw and rotate model using 3D scatter in `matplotlib` with `%matplotlib notebook`.\n", "\n", "\n", "* Download the 3D model of the bunny [here](https://drive.google.com/open?id=1H_pqd2fInkkCKF5R9voVFJM53233G-Ax).\n", "\n", "\n", "* You can read it and display with `mesh = trimesh.load()` and `mesh.show()` respectively. Coordinates of the vertices and faces can be accessed with `mesh.vertices` and `mesh.faces` respectively. Note that not every vertex is a part of a face, you should filter them out.\n", "\n", "\n", "* Create an adjacency matrix of the graph corresponding to the mesh of the model (you can use `networkx` library for this). Verify that the number of connected components of the graph is exactly $1$.\n", "\n", "\n", "* Find the \"bunny\" partitioning for `num_iter_fix=1, num_iter_adapt=2` and `num_iter_fix=2, num_iter_adapt=10`. Plot the convergence rate for the second case and discuss it. Draw 3D model with coloring according to the partitions. With `trimesh` you can assign vertex colors as `mesh.visual.vertex_colors[] = (R,G,B,A)`." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem 3 (30 pts)\n", "\n", "## Sparse matrices and linear solvers" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consider a 2D convection-diffusion equation in $\\Omega = [0,1]^2$\n", "$$\n", "-\\frac{\\partial^2 u}{\\partial x^2} - \\frac{\\partial^2 u}{\\partial y^2} + \\frac{\\partial u}{\\partial x} + \\frac{\\partial u}{\\partial y} = f(x,y), \\quad (x,y)\\in \\Omega\n", "$$ \n", "with zero Dirichlet boundary conditions\n", "$$\n", "u_{\\partial \\Omega} = 0,\n", "$$\n", "with known function $f(x,y)$ and unknown $u(x,y)$.\n", "\n", "To find solution of this equation we will use the **finite difference method**. \n", "Standard second order finite difference discretization on a uniform grid $(x_i, y_j) = (ih, jh)$, $i,j = 0,\\dots, N$, $h = \\frac{1}{N}$ leads to the following system of equations:\n", "$$\n", "\\begin{split}\n", " &-\\frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{h^2} - \\frac{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{h^2} + \\frac{u_{i+1, j} - u_{i,j}}{h} + \\frac{u_{i, j+1} - u_{i,j}}{h} = f(ih, jh) \\\\\n", " &u_{0,j} = u_{i,0} = u_{N,j} = u_{i,N} = 0, \\quad i,j = 0,\\dots,N,\n", "\\end{split}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Write the system above as a matrix equation $BU_h + U_h C = F_h$ with matrices $U_h = \\begin{bmatrix}u_{1,1} & \\dots & u_{1,N-1} \\\\ \\vdots & \\ddots & \\vdots \\\\ u_{N-1,1} & \\dots & u_{N-1,N-1} \\end{bmatrix}$, $F_h = \\begin{bmatrix}f_{1,1} & \\dots & f_{1,N-1} \\\\ \\vdots & \\ddots & \\vdots \\\\ f_{N-1,1} & \\dots & f_{N-1,N-1} \\end{bmatrix}$. What are matrices $B$ and $C$?\n", "\n", "* Using Kronecker product properties rewrite **analytically** $ BU_h + U_h C = F_h$ as $A_h \\mathrm{vec}(U_h) = \\mathrm{vec}(F_h)$, where $\\mathrm{vec}(\\cdot)$ is a columnwise reshape.\n", "What is matrix $A_h$?\n", "\n", "* Choose $f(x,y) = 1$. \n", " Solve the system with this $f$ using the `scipy.sparse.linalg.spsolve` which is direct sparse solver. \n", " Use ```pandas``` library and print table that contains $N$ and time of solving for $N=127, 255, 511$. \n", " Matrices $B, C$ and $A_h$ should be assembled in the `CSR` format using functions from the `scipy.sparse` package (functions `scipy.sparse.kron` and `scipy.sparse.spdiags` will be helpful). Do not use full matrices! Use only sparse arithmetics. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# INPUT - dimension of the grid along the one side of the square\n", "# OUTPUT - matrix Ah for 2D convection-diffusion problem in CSR sparse format \n", "def build_matrix(N): # 5 pts\n", " # Your code is here\n", " return A" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Fix initial random guess $x_0$, maximal number of iterations and required tolerance. These parameters will be the same for all iterative methods that you will test\n", "* Run `cg`, `minres`, `GMRES`, `BicgStab` from `scipy.sparse.linalg` package with generated matrix for 2D convection-diffusion equation, right-hand side $f(x, y) = 1$ and above described parameters.\n", "* What are the iterative methods diverge and why?\n", "* What are the iterative methods of choice? Explain why. \n", "\n", "* Plot the relative residual norm $\\frac{\\|Ax_k - f\\|_2}{\\|f\\|_2}$ w.r.t number of iteration and w.r.t elapsed time. Assume that every iteration lasts the same number of seconds\n", "* Run the methods of choice with and without ILU0 preconditioner (use it from `scipy.sparse` package) and Block Jacobi preconditioner for $N=256$. You should implement Block Jacobi preconditioner in the most efficient manner, in particular the most complex operation is solve linear systems with $N \\times N$ matrix and multiple right-hand sides. Plot relative error w.r.t iteration number for both cases and for both preconditioners on one plot. Don't forget add legend. Also plot relative error w.r.t time in the similar way. Assume that every iteration lasts the same number of seconds" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Problem 4 (30 pts)\n", "\n", "## Again QR, oh no!\n", "\n", "\n", "* In this problem you asked to find the convolution of the $n\\times n$, $n=512$ QR(-code) with the following filter\n", "$$\n", "T_{i_1j_1,i_2j_2} \\equiv T_{i_1-j_1,i_2-j_2} = \\frac{\\alpha}{\\pi} e^{-\\alpha \\left[(i_1 - j_1)^2 + (i_2 - j_2)^2 \\right]}, \\quad i_1,j_1, i_2, j_2 = 1,\\dots, n, \\quad \\alpha = \\frac{1}{100}\n", "$$\n", "using FFT. \n", " * Write function `matvec` that produces multiplication of $T$ by a given vector $x$. \n", " * Use `scipy.sparse.linalg.LinearOperator` to create an object that has attribute `.dot()` (this object will be further used in the iterative process). Note that `.dot()` input and output must be 1D vectors, so do not forget to use reshape.\n", " * Convolve `qrcode.jpg` with $T$ for $\\alpha = \\frac{1}{100}$.Plot the result as an image. \n", " * What is the complexity of this operation? \n", " *\n", "**Note:** You can use standart FFT, from e.q. `numpy.fft`\n", "\n", "\n", "* Run an appropriate Krylov method(s) with the obtained Linear Operator and try to reconstruct back QR-code using the right-hand side from the first bullet (smoothed QR-code). \n", "* On one figure plot norm of residual with respect to the number of iterations for $\\alpha = \\{\\frac{1}{50},~~\\frac{1}{100},~~\\frac{1}{200}\\}$ and corresponding right hand side. Comment on the results.\n", "\n", "\n", "* **Bonus:** Find image `noisy_extra.jpg` in the attached to this problem set file. This image is the QR-code with the beautiful message, convolved with $T$, parameter $\\alpha=\\frac{1}{100}$ plus small random noise:\n", "\n", " * Your goal now is to get the link from the image. Since the problem is very ill-posed you need some regularization, denoising or anything else." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#Function that will be provided into scipy.sparse.linalg.LinearOperator\n", "#INPUT: 1D array %vec and \n", "# alpha with default value==1/100.\n", "#OUTPUT: 1D array\n", "#Hint: you can vary alpha using lambda function as argument of linalg.LinearOperator:\n", "# e.q. lambda x: matvec(x,1/50.)\n", "def matvec(vec, alpha=1/100.): # 15 pts\n", " return out" ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3", "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.6.3" } }, "nbformat": 4, "nbformat_minor": 2 }