{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# A much too brief tour of numerical linear algebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We want to compute an eigenvalue $\\lambda$ of a square matrix $A$ and its eigenvector $x$, i.e.\n",
    "$$\n",
    "Ax=\\lambda x\n",
    "$$\n",
    "using the power method (Google's page rank used this)\n",
    "1. Choose $q^0$\n",
    "2. Compute $q^{(k+1)}=\\frac{Aq^{(k)}}{\\Vert Aq^{(k)}\\Vert}$, **Repeat this step until convergence**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Computing an eigenvalue and eigenvector of a matrix\n",
    "We assume that $A$ is diagonalizable,\n",
    "\\begin{align*}X^{-1} AX\n",
    "  =\\begin{bmatrix}\n",
    "    \\lambda_1&~&~\\\\~&\\ddots&~\\\\~&~&\\lambda_n\n",
    "  \\end{bmatrix}\n",
    "\\end{align*}\n",
    "with $X=\\left[ x_1,\\dots, x_n\\right], \\; \\|x_j\\|_2=1,\\; |\\lambda_{1}| > |\\lambda_2|\n",
    "\\geq \\dots \\geq |\\lambda_n|$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Under the above assumptions let\n",
    "  \\begin{align*}\n",
    "    A x_j =\\lambda_j x_j, ~\\text{und}~ q^{(0)}=\\sum\\limits_{j=1}^n \\alpha_j x_j \n",
    "  \\end{align*}\n",
    "  with $\\alpha_1\\neq 0$. Then\n",
    "  \\begin{align*}\n",
    "    \\lim\\limits_{k \\rightarrow \\infty} \\lambda^{(k)}&= \\lambda_1,~\\lim\\limits_{k\n",
    "      \\rightarrow \\infty} q^{(k)} = \\mu \\cdot x_1,\\\\\n",
    "    |\\mu|&=1\\quad(\\Rightarrow \\mu = e^{i\\omega}, \\omega \\in\\mathbb{R})\n",
    "  \\end{align*}\n",
    "  with convergence rate $\\sim |\\lambda_2| / |\\lambda_1| < 1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We start with\n",
    "$$\n",
    "q^{(k)}=\\gamma_k A q^{(k-1)} =(\\gamma_k\\gamma_{k-1} A^{2}q^{(k-2)}) =\\dots\n",
    "    =\\underbrace{\\gamma_k \\dots \\gamma_1}_{=:\\beta_k} A^kq^{(0)}\n",
    "$$\n",
    "and have\n",
    "$$  \n",
    "  \\gamma_k=1 / \\|A q^{(k-1)}\\|_2\\\\\n",
    "$$\n",
    "$$\n",
    " q^{(k)}=\\beta_k\\;A^k q^{(0)},~ \\text{also}~ q^{(k)} \\in \\mathrm{span}{A^k q^{(0)}}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "It holds that\n",
    "  \\begin{align*}\n",
    "    \\frac{A^k q^{(0)}}{\\Vert A^k q^{(0)}\\Vert}&=\\frac{A^k \\sum\\limits_{j=1}^n \\alpha_j x_j}{\\Vert A^k q^{(0)}\\Vert}\\\\\n",
    "    &=\\frac{A^{(k-1)} \\sum\\limits_{j=1}^n \\alpha_j \\underbrace{A\\; x_j}_{\\lambda_j x_j}}{\\Vert A^k q^{(0)}\\Vert}\\\\\n",
    "    &=\\frac{\\alpha_1 \\lambda_1^k \\left(x_1 + \\sum\\limits_{j=2}^n \\frac{\\alpha_j}{\\alpha_1}\n",
    "      \\left(\\frac{\\lambda_j}{\\lambda_1} \\right)^k \\; x_j\\right)}{\\vert\\alpha_1\\vert \\vert\\lambda_1\\vert^k \\Vert x_1 + \\sum\\limits_{j=2}^n \\frac{\\alpha_j}{\\alpha_1}\n",
    "      \\left(\\frac{\\lambda_j}{\\lambda_1} \\right)^k \\; x_j\\Vert}\n",
    "  \\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "  and so the limit is \n",
    "  $$\n",
    "  \\frac{\\alpha_1 \\lambda_1^k x_1}{\\vert\\alpha_1\\vert \\vert\\lambda_1\\vert^k}\n",
    "  $$\n",
    "  $\\Rightarrow \\mathrm{dist}(\\mathrm{span}{q^{(k)}},~\\mathrm{span}{x_1}) \\rightarrow 0$,\n",
    "  and the convergence speed is then given by\n",
    "  $\\frac{|\\lambda_2|}{|\\lambda_1|}. \\quad $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's see whether I can implement the power method myself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from scipy import linalg as la\n",
    "from scipy import optimize\n",
    "import scipy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's learn about the QR factorisation where we can factor a matrix\n",
    "$$\n",
    "A\\in\\mathbb{R}^{m,n}\n",
    "$$\n",
    "is decomposed as\n",
    "$$\n",
    "A=QR\n",
    "$$\n",
    "with $Q\\in\\mathbb{R}^{m,m}$ orthogonal\n",
    "and \n",
    "$$\n",
    "R=\n",
    "\\begin{bmatrix}\n",
    "\\hat{R}\\\\\n",
    "0\\\\\n",
    "\\end{bmatrix}\\in\\mathbb{R}^{m,n}\n",
    "$$\n",
    "and \n",
    "$\\hat{R}\\in\\mathbb{R}^{n,n}$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "la.qr?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "n=10\n",
    "A=np.random.rand(n,n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.27086569 0.80317013 0.79002127 0.18604938 0.5868645  0.57189462\n",
      "  0.76964117 0.28441925 0.42529308 0.79435346]\n",
      " [0.47789391 0.98695112 0.31192841 0.73986469 0.01808935 0.07802996\n",
      "  0.95102272 0.35596123 0.74238573 0.10866131]\n",
      " [0.65065064 0.3669956  0.31034799 0.43518505 0.86478858 0.34963775\n",
      "  0.36001982 0.4152881  0.54611148 0.8551706 ]\n",
      " [0.05559374 0.87812647 0.07804075 0.38475161 0.92414789 0.41092123\n",
      "  0.96496381 0.1492922  0.5889473  0.19987846]\n",
      " [0.52380139 0.58418249 0.73158582 0.98092487 0.18787021 0.41200372\n",
      "  0.30963949 0.54824776 0.87469683 0.19889076]\n",
      " [0.1093706  0.32248406 0.44152885 0.97185763 0.19552999 0.12977521\n",
      "  0.76455776 0.75936689 0.96977646 0.93969089]\n",
      " [0.75979629 0.46713386 0.37077584 0.95642245 0.57309286 0.51880711\n",
      "  0.70925084 0.84958429 0.58798348 0.80735443]\n",
      " [0.38734541 0.10653562 0.49686729 0.9655784  0.22938752 0.87835686\n",
      "  0.9001694  0.56243919 0.93636609 0.16723133]\n",
      " [0.08101855 0.22833441 0.40125029 0.96972837 0.9733203  0.55680207\n",
      "  0.13471853 0.41314343 0.61112646 0.59818568]\n",
      " [0.89057594 0.02766626 0.52509287 0.34391331 0.93811894 0.27376685\n",
      "  0.89857892 0.67648181 0.87030987 0.55132657]]\n",
      "[[0.27086569 0.80317013 0.79002127 0.18604938 0.5868645  0.57189462\n",
      "  0.76964117 0.28441925 0.42529308 0.79435346]\n",
      " [0.47789391 0.98695112 0.31192841 0.73986469 0.01808935 0.07802996\n",
      "  0.95102272 0.35596123 0.74238573 0.10866131]\n",
      " [0.65065064 0.3669956  0.31034799 0.43518505 0.86478858 0.34963775\n",
      "  0.36001982 0.4152881  0.54611148 0.8551706 ]\n",
      " [0.05559374 0.87812647 0.07804075 0.38475161 0.92414789 0.41092123\n",
      "  0.96496381 0.1492922  0.5889473  0.19987846]\n",
      " [0.52380139 0.58418249 0.73158582 0.98092487 0.18787021 0.41200372\n",
      "  0.30963949 0.54824776 0.87469683 0.19889076]\n",
      " [0.1093706  0.32248406 0.44152885 0.97185763 0.19552999 0.12977521\n",
      "  0.76455776 0.75936689 0.96977646 0.93969089]\n",
      " [0.75979629 0.46713386 0.37077584 0.95642245 0.57309286 0.51880711\n",
      "  0.70925084 0.84958429 0.58798348 0.80735443]\n",
      " [0.38734541 0.10653562 0.49686729 0.9655784  0.22938752 0.87835686\n",
      "  0.9001694  0.56243919 0.93636609 0.16723133]\n",
      " [0.08101855 0.22833441 0.40125029 0.96972837 0.9733203  0.55680207\n",
      "  0.13471853 0.41314343 0.61112646 0.59818568]\n",
      " [0.89057594 0.02766626 0.52509287 0.34391331 0.93811894 0.27376685\n",
      "  0.89857892 0.67648181 0.87030987 0.55132657]]\n"
     ]
    }
   ],
   "source": [
    "q,r = la.qr(A)\n",
    "print(A)\n",
    "print(q@r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
       "       [ 0.,  2.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  3.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  4.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  0.,  5.,  0.,  0.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  0.,  0.,  6.,  0.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  0.,  0.,  0.,  7.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  8.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  9.,  0.],\n",
       "       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., 10.]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n=10\n",
    "D=np.diag(np.linspace(1,n,n))\n",
    "D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's create a matrix with given eigenvalues!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "AA=np.dot(np.transpose(q),np.dot(D,q)) #a similarity transformation AA=Q^TDQ with D containing the eigenvalues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 6.48734266 -1.91569246 -0.33522128 -0.23752101 -0.80076558  0.09716313\n",
      "   1.291004    0.0298566  -0.45188107 -1.62860994]\n",
      " [-1.91569246  4.24668069  0.11604303 -1.64938555 -0.01770765 -0.01462223\n",
      "  -1.36395159 -0.62774587  0.74901977  0.28935442]\n",
      " [-0.33522128  0.11604303  4.60208768  2.1615595   0.39038219 -0.28612837\n",
      "  -0.1573882  -1.71442713 -0.03380141 -0.44568929]\n",
      " [-0.23752101 -1.64938555  2.1615595   5.13181275  0.60934557  1.23437768\n",
      "   0.78204937  0.84974057  1.00294252 -1.01865992]\n",
      " [-0.80076558 -0.01770765  0.39038219  0.60934557  5.59488307 -1.24148594\n",
      "   0.16402429 -1.40917157  1.8936234   0.21826504]\n",
      " [ 0.09716313 -0.01462223 -0.28612837  1.23437768 -1.24148594  7.09066893\n",
      "  -0.02175411 -0.34353213 -0.2727222  -0.55818141]\n",
      " [ 1.291004   -1.36395159 -0.1573882   0.78204937  0.16402429 -0.02175411\n",
      "   6.48378605  0.27067475  1.02298559  0.12844229]\n",
      " [ 0.0298566  -0.62774587 -1.71442713  0.84974057 -1.40917157 -0.34353213\n",
      "   0.27067475  5.98592775  0.07005447 -0.08681968]\n",
      " [-0.45188107  0.74901977 -0.03380141  1.00294252  1.8936234  -0.2727222\n",
      "   1.02298559  0.07005447  5.01893767  0.19085268]\n",
      " [-1.62860994  0.28935442 -0.44568929 -1.01865992  0.21826504 -0.55818141\n",
      "   0.12844229 -0.08681968  0.19085268  4.35787276]]\n"
     ]
    }
   ],
   "source": [
    "print(AA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.51631768],\n",
       "       [0.71709195],\n",
       "       [0.27783863],\n",
       "       [0.4420662 ],\n",
       "       [0.45700284],\n",
       "       [0.40666308],\n",
       "       [0.25739401],\n",
       "       [0.93299313],\n",
       "       [0.17620844],\n",
       "       [0.46853357]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xold=np.random.rand(10,1)\n",
    "xold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[9.99998333]]\n"
     ]
    }
   ],
   "source": [
    "itmax = 50\n",
    "xold=np.random.rand(10,1)\n",
    "# main iteration loop for power method\n",
    "for n in range(itmax):\n",
    "    y=np.dot(AA,xold)\n",
    "    xnew=y/la.norm(y)\n",
    "    lmbd=(np.dot(np.transpose(xnew),np.dot(AA,xnew))) # What does this look like of xnew=x1?\n",
    "    xold=xnew\n",
    "    # Here comes your stopping criteria\n",
    "    \n",
    "print(lmbd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "la.norm?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Solving a linear system\n",
    "\n",
    "We are now interested in solving a linear system of equations, which could be obtained from the finite element discretization of the partial differential equation \n",
    "$$\n",
    "-\\Delta u=f\\quad\\mathrm{on}\\quad\\Omega\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "with appropriate boundary conditions. And the result is\n",
    "$$\n",
    "\\frac{1}{h^2}\n",
    "\\left[\n",
    "\\begin{array}{cccc}\n",
    "2&-1&&\\\\\n",
    "-1&2&\\ddots&\\\\\n",
    "&\\ddots&\\ddots&-1\\\\\n",
    "&&-1&2\\\\\n",
    "\\end{array}\n",
    "\\right]\n",
    "\\left[\n",
    "\\begin{array}{c}\n",
    "u_1\\\\\n",
    "u_2\\\\\n",
    "\\vdots\\\\\n",
    "u_n\\\\\n",
    "\\end{array}\n",
    "\\right]\n",
    "=\n",
    "\\left[\n",
    "\\begin{array}{c}\n",
    "f_1\\\\\n",
    "f_2\\\\\n",
    "\\vdots\\\\\n",
    "f_n\\\\\n",
    "\\end{array}\n",
    "\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "with the system matrix\n",
    "$$\n",
    "A=\n",
    "\\frac{1}{h^2}\\left[\n",
    "\\begin{array}{cccc}\n",
    "2&-1&&\\\\\n",
    "-1&2&\\ddots&\\\\\n",
    "&\\ddots&\\ddots&-1\\\\\n",
    "&&-1&2\\\\\n",
    "\\end{array}\n",
    "\\right],\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "right hand side vector\n",
    "$$\n",
    "f=\n",
    "\\left[\n",
    "\\begin{array}{c}\n",
    "f_1\\\\\n",
    "f_2\\\\\n",
    "\\vdots\\\\\n",
    "f_n\\\\\n",
    "\\end{array}\n",
    "\\right]\n",
    "$$\n",
    "and solution vector\n",
    "$$\n",
    "u=\n",
    "\\left[\n",
    "\\begin{array}{c}\n",
    "u_1\\\\\n",
    "u_2\\\\\n",
    "\\vdots\\\\\n",
    "u_n\\\\\n",
    "\\end{array}\n",
    "\\right].\n",
    "$$\n",
    "So we now want so solve the system \n",
    "$$\n",
    "Au=f.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.sparse as spsp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "N = 6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "diagonals = np.zeros((3, N))   # 3 diagonals\n",
    "diagonals[0,:] = -1\n",
    "diagonals[1,:] = 2\n",
    "diagonals[2,:] = -1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "A = spsp.spdiags(diagonals, [-1,0,1], N, N, format='csc')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  (1, 0)\t-1.0\n",
      "  (0, 0)\t2.0\n",
      "  (2, 1)\t-1.0\n",
      "  (1, 1)\t2.0\n",
      "  (0, 1)\t-1.0\n",
      "  (3, 2)\t-1.0\n",
      "  (2, 2)\t2.0\n",
      "  (1, 2)\t-1.0\n",
      "  (4, 3)\t-1.0\n",
      "  (3, 3)\t2.0\n",
      "  (2, 3)\t-1.0\n",
      "  (5, 4)\t-1.0\n",
      "  (4, 4)\t2.0\n",
      "  (3, 4)\t-1.0\n",
      "  (5, 5)\t2.0\n",
      "  (4, 5)\t-1.0\n"
     ]
    }
   ],
   "source": [
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-1.  -0.6 -0.2  0.2  0.6  1. ]\n"
     ]
    }
   ],
   "source": [
    "x = np.linspace(-1, 1, N) \n",
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-1.40000000e+00  0.00000000e+00 -1.11022302e-16  2.22044605e-16\n",
      "  0.00000000e+00  1.40000000e+00]\n"
     ]
    }
   ],
   "source": [
    "b = A.dot(x)\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "expected matrix",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-19-c3fc74ac55c0>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;31m# Perform Gaussian elimination A=PLU\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mP\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mL\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mU\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mla\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlu\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mA\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/miniconda3/envs/scp/lib/python3.7/site-packages/scipy/linalg/decomp_lu.py\u001b[0m in \u001b[0;36mlu\u001b[0;34m(a, permute_l, overwrite_a, check_finite)\u001b[0m\n\u001b[1;32m    212\u001b[0m         \u001b[0ma1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0masarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    213\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 214\u001b[0;31m         \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'expected matrix'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    215\u001b[0m     \u001b[0moverwrite_a\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0moverwrite_a\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0m_datacopied\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    216\u001b[0m     \u001b[0mflu\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_flinalg_funcs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'lu'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0ma1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: expected matrix"
     ]
    }
   ],
   "source": [
    "# Perform Gaussian elimination A=PLU\n",
    "P, L, U = la.lu(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "la.lu?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "So that did not work with our sparse matrix format. What else can we do?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "factor=spsp.linalg.splu(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's see what we have got"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<SuperLU object at 0x7f4d9cd95ed0>\n"
     ]
    }
   ],
   "source": [
    "print(factor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-1.  -0.6 -0.2  0.2  0.6  1. ]\n"
     ]
    }
   ],
   "source": [
    "x=factor.solve(b)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-1.  -0.6 -0.2  0.2  0.6  1. ]\n"
     ]
    }
   ],
   "source": [
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.925231146709438e-17"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scipy.linalg.norm(A*x-b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "This previous working is based on the factorization of $A$ into two triangular factors and tybpically a permutation matrix $P$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "spsp.linalg?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}