{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 9: Matrix decompositions and how we compute them"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Syllabus\n",
    "**Week 1:** Matrices, vectors, matrix/vector norms, scalar products & unitary matrices  \n",
    "**Week 2:** TAs-week (Strassen, FFT, a bit of SVD)  \n",
    "**Week 3:** Matrix ranks, singular value decomposition, linear systems, eigenvalues  \n",
    "**Week 4:** Matrix decompositions: QR, LU, SVD + test + structured matrices start"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Recap of the previous lecture\n",
    "- Eigenvectors and eigenvalues\n",
    "- Characteristic polynomial and why it is a bad idea\n",
    "- Power method: $x := Ax$\n",
    "- Gershgorin theorem\n",
    "- Schur theorem: $A = U T U^*$ \n",
    "- Normal matrices: $A^* A = A A^*$\n",
    "- A brief whiteboard illustration of how to compute the Schur form: QR-algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Today lecture\n",
    "Today we will talk about **matrix factorizations**\n",
    "\n",
    "- LU decomposition and Gaussian elimination\n",
    "- QR decomposition and Gram-Schmidt algorithm\n",
    "- Schur decomposition and QR-algorithm\n",
    "\n",
    "We already introduced it some time ago, but now we are going to do it in more details.  \n",
    "\n",
    "These are the basic **matrix factorizations** in numerical linear algebra."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## General concept of matrix factorization\n",
    "\n",
    "In numerical linear algebra we need to solve different tasks, for example:\n",
    "\n",
    "- Solve linear systems $Ax = f$\n",
    "- Compute eigenvalues / eigenvectors\n",
    "- Compute singular values / singular vectors\n",
    "- Compute inverses, even sometimes determinants \n",
    "- Compute **matrix functions** like $\\exp(A), \\cos(A)$.\n",
    "\n",
    "In order to do this, we replace the matrix as a sum and/or product of matrices with **simpler structure**,  \n",
    "such that we can solve abovementioned tasks faster / more stable form.\n",
    "\n",
    "What is a **simpler structure**?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## What is a simpler structure\n",
    "We already encountered several classes of matrices **with structure**. \n",
    "\n",
    "In dense matrix business, the most important classes are\n",
    "\n",
    "**unitary matrices** ($U^* U = I$, why?), **upper/lower triangular matrices **, **diagonal matrices**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Other classes of matrices\n",
    "\n",
    "For **sparse matrices** the sparse constraints are often included in the factorizations.  \n",
    "\n",
    "For **Toeplitz matrices** an important class of matrices is the matrices with **low displacement rank**, \n",
    "\n",
    "which is based on the low-rank matrices. \n",
    "\n",
    "The class of **low-rank matrices** and **block low-rank matrices** is everywhere."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Plan\n",
    "The plan for todays lecture is to discuss the decompositions one-by-one and point out:\n",
    "- How to compute a particular decomposition\n",
    "- When the decomposition exists \n",
    "- What is done in \"real life\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Decompositions we want to discuss today\n",
    "- LU factorization\n",
    "- Positive definite matrices and Cholesky factorization\n",
    "- QR-decomposition and Gram-Schmidt algorithm\n",
    "- 1 slide about the SVD (more tomorrow)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## What is the LU factorization\n",
    "\n",
    "LU factorization is the representation of a given matrix $A$ as a product  \n",
    "of a **lower triangular matrix** $L$ and an **upper triangular** matrix $U$:\n",
    "\n",
    "$$A = LU$$\n",
    "\n",
    "This factorization is **non-unique**, so it is typical to require that the matrix $L$ has ones on the diagonal.\n",
    "\n",
    "Does it always exist?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Main goal** of the LU decomposition is to solve linear systems, because\n",
    "\n",
    "$$\n",
    "    A^{-1} f = (L U)^{-1} f = U^{-1} L^{-1} f, \n",
    "$$\n",
    "\n",
    "and this reduces to the solution of two linear systems \n",
    "$$\n",
    "     L y = f, \n",
    "$$\n",
    "and \n",
    "$$\n",
    "   U x = y.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Solving linear systems with **triangular matrices** is **easy**: you just solve the equation with one unknown, put into the second equation, and so on.  \n",
    "\n",
    "The computational cost of such algorithm is $\\mathcal{O}(n^2)$ arithmetic operations, once $LU$ factorization is known. \n",
    "\n",
    "Note that the **inverse of triangular matrix** is a triangular matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "⎡ 1    0    0   0   0⎤\n",
      "⎢                    ⎥\n",
      "⎢1/3   1    0   0   0⎥\n",
      "⎢                    ⎥\n",
      "⎢1/5  6/7   1   0   0⎥\n",
      "⎢                    ⎥\n",
      "⎢          15        ⎥\n",
      "⎢1/7  5/7  ──   1   0⎥\n",
      "⎢          11        ⎥\n",
      "⎢                    ⎥\n",
      "⎢     20   210  28   ⎥\n",
      "⎢1/9  ──   ───  ──  1⎥\n",
      "⎣     33   143  15   ⎦\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import sympy\n",
    "from sympy.matrices import Matrix\n",
    "import IPython\n",
    "from sympy.interactive.printing import init_printing\n",
    "init_printing(use_latex=True)\n",
    "n = 5\n",
    "w = Matrix(n, n, lambda i, j: 1/(i + j +  sympy.Integer(1)/2))\n",
    "L, U, tmp = w.LUdecomposition()\n",
    "#Generate the final expression\n",
    "#fn = u'%s \\\\times %s = %s' % (sympy.latex(L), sympy.latex(U), sympy.latex(w))\n",
    "#L * U - w\n",
    "L"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The classical LU-decomposition algorithm\n",
    "The simplest way to implement LU-decomposition is the following 3-cycle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from sympy.interactive.printing import init_printing\n",
    "init_printing(use_latex=True)\n",
    "L = Matrix(n, n, lambda i, j: 0)\n",
    "U = Matrix(n, n, lambda i, j: 0)\n",
    "a = w.copy() #Our matrix\n",
    "for k in xrange(n): #Eliminate one row\n",
    "    L[k, k] = 1\n",
    "    for i in xrange(k+1, n):\n",
    "        L[i, k] = a[i, k] / a[k, k]\n",
    "        for j in xrange(k+1, n):\n",
    "            a[i, j] = a[i, j] - L[i, k] * a[k, j]\n",
    "    for j in xrange(k, n):\n",
    "        U[k, j] = a[k, j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Existence of the LU-decomposition\n",
    "The LU-decomposition algorithm does not fail  \n",
    "\n",
    "if **we do not divide by zero** at every step.\n",
    "\n",
    "When it is so, for which class of matrices?\n",
    "\n",
    "It is true for **strictly regular matrices**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Strictly regular matrices and LU decomposition \n",
    "A matrix $A$ is called strictly regular, if all of its **principal minors**  \n",
    "(i.e, submatrices in the first $k$ rows and $k$ columns) are non-singular. \n",
    "\n",
    "In this case, there always exists a LU-decomposition.  The reverse is also true (check!)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Proof**. If there is a LU-decomposition, then the matrix is strictly regular \n",
    "\n",
    "This follows from the fact that to get a minor you multiply a corresponding submatrix of $L$ by a corresponding submatrix of $U$, and they are non-singular (since non-singularity of triangular matrices is equivalent to the fact that their diagonal elements are not equal to zero).\n",
    "\n",
    "The other way can be proven by **induction**. Suppose that we know that for all $(n-1) \\times (n-1)$ matrices will non-singular minors LU-decomposition exists."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Then, consider the block form\n",
    "$$\n",
    "    A = \\begin{bmatrix}\n",
    "          a & c^{\\top} \\\\\n",
    "          b & D\n",
    "    \\end{bmatrix},\n",
    "$$\n",
    "where $D$ is $(n-1) \\times (n-1)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Then we do \"block elimination\": \n",
    "\n",
    "$$\n",
    "     \\begin{bmatrix}\n",
    "     1 & 0 \\\\\n",
    "     -z & I\n",
    "     \\end{bmatrix}\n",
    "     \\begin{bmatrix}\n",
    "     a & c^{\\top} \\\\\n",
    "     b & D \n",
    "     \\end{bmatrix}=\n",
    "     \\begin{bmatrix}\n",
    "     a & c^{\\top} \\\\\n",
    "     0 & A_1\n",
    "     \\end{bmatrix},\n",
    "$$\n",
    "\n",
    "where $z = \\frac{b}{a}, \\quad A_1 = D - \\frac{1}{a} b c^{\\top}$.  \n",
    "We can show that $A_1$ is also strictly regular, thus it has (by induction) the LU-decomposition:  \n",
    "$$\n",
    "   A_1 = L_1 U_1, \n",
    "$$\n",
    "And that gives the $LU$ decomposition of the original matrix.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## When LU fails\n",
    "What happens, if the matrix is not strictly regular (or the **pivots** in the Gaussian elimination are really small?). \n",
    "\n",
    "There is classical $2 \\times 2$ example of a matrix with a bad LU decomposition.  \n",
    "\n",
    "The matrix we look at is  \n",
    "$$\n",
    "    A = \\begin{pmatrix}\n",
    "    \\varepsilon & 1 \\\\\n",
    "    1 & 1 \n",
    "    \\end{pmatrix}\n",
    "$$\n",
    "\n",
    "If $\\varepsilon$ is sufficiently small, we **might** fail.\n",
    "\n",
    "Let us do some demo here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "L * U - A: [[ 0.  0.]\n",
      " [ 0.  0.]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[  1.00000000e+00,   0.00000000e+00],\n",
       "       [  8.92857143e+15,   1.00000000e+00]])"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "eps = 1.12e-16\n",
    "a = [[eps, 1],[1.0,  1]]\n",
    "a = np.array(a)\n",
    "a0 = a.copy()\n",
    "n = a.shape[0]\n",
    "L = np.zeros((n, n))\n",
    "U = np.zeros((n, n))\n",
    "for k in xrange(n): #Eliminate one row   \n",
    "    L[k, k] = 1\n",
    "    for i in xrange(k+1, n):\n",
    "        L[i, k] = a[i, k] / a[k, k]\n",
    "        for j in xrange(k+1, n):\n",
    "            a[i, j] = a[i, j] - L[i, k] * a[k, j]\n",
    "    for j in xrange(k, n):\n",
    "        U[k, j] = a[k, j]\n",
    "\n",
    "print 'L * U - A:', np.dot(L, U) - a0\n",
    "L"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The concept of pivoting\n",
    "We can do pivoting, i.e. permute rows and columns to maximize $A_{kk}$ that we divide over.  \n",
    "The simplest but effective strategy is the **row pivoting**: at each step, select the index that is maximal in modulus, and put it onto the diagonal. \n",
    "\n",
    "It gives us the decomposition \n",
    "\n",
    "$$A = P L U,$$\n",
    "where $P$ is a **permutation matrix**.\n",
    "\n",
    "\n",
    "What makes **row pivoting** good?  \n",
    "\n",
    "It is made good by the fact that  \n",
    "\n",
    "$$\n",
    "   | L_{ij}|<1,\n",
    "$$\n",
    "but the elements of $U$ can grow, up to $2^n$! (in practice, this is very rarely encountered).\n",
    "\n",
    "Can you come up with a matrix where the elements of $U$ grow as much as possible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Positive definite matrices\n",
    "\n",
    "Strictly regular matrices have LU-decomposition. \n",
    "\n",
    "An important **subclass** of strictly regular matrices is the class of **Hermitian positive definite matrices**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Positive definite matrices(2)\n",
    "**Definition:** A matrix is called <font color='red'> positive definite </font> if for any $x: \\Vert x \\Vert \\ne 0$ we have\n",
    "\n",
    "$$(Ax, x) > 0.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Symmetric positive definite matrices\n",
    "\n",
    "**Statement:** A Hermitian positive definite matrix $A$ is strictly regular and has **Cholesky factorization** of the form \n",
    "\n",
    "$$A = RR^*,$$\n",
    "\n",
    "where $R$ is upper triangular.\n",
    "\n",
    "Let us try to prove this fact (on the whiteboard).\n",
    "\n",
    "The Cholesky factorization is **always stable**.\n",
    "\n",
    "It is sometimes referred to as \"square root\" of the matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Summary on the LU for dense matrices\n",
    "- Principal submatrices non-singular -- there exists LU\n",
    "- Pivoting is needed to avoid error growth\n",
    "- Complexity is $\\mathcal{O}(n^3)$ for the factorization step (can be speeded by using Strassen ideas) and $\\mathcal{O}(n^2)$ for a solve\n",
    "- For SPD matrices Cholesky factorization is the method of choice (stability, twice less memory to store)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR-decomposition\n",
    "The next decomposition: **QR** decomposition. Again from the name it is clear that a matrix is represented as a product  \n",
    "$$\n",
    "    A = Q R, \n",
    "$$\n",
    "where $Q$ is an **orthogonal (unitary)** matrix and $R$ is **upper triangular**.  \n",
    "\n",
    "The matrix sizes: $Q$ is $n \\times m$, $R$ is $m \\times m$.\n",
    "\n",
    "QR-factorization is defined for **any rectangular matrix**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR-decomposition: applications\n",
    "\n",
    "This algorithm plays a crucial role in many problems:\n",
    "- Computing orthogonal bases in a linear space\n",
    "- Used in the preprocessing step for the SVD\n",
    "- QR-algorithm for the computation of eigenvectors and eigenvalues (one of the 10 most important algorithms of the 20th century) is    based on the QR-decomposition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Theorem \n",
    "Every rectangular $n \\times m$ matrix, $n \\geq m$ matrix has a QR-decomposition. \n",
    "There are several ways to prove it and compute it.\n",
    "\n",
    "- (theoretical) Using the Gram matrices and Cholesky factorization\n",
    "- (geometrical) Using the Gram-Schmidt orthogonalization\n",
    "- (practical) Using Householder/Givens transformations \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Mathematical way\n",
    "If we have the representation of the form\n",
    "$$A = QR,$$\n",
    "then $A^* A = R^* R$, the matrix $A^* A$ is called **Gram matrix**, and its elements are scalar products of the columns of $A$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Math way: full column rank\n",
    "\n",
    "Assume that $A$ has **full column rank**. Then, \n",
    "\n",
    "It is simple to show that $A^* A$ is positive definite:\n",
    "\n",
    "$$\n",
    "   (A^* A y, y) = (Ay, Ay)^2 = \\Vert Ay \\Vert  \\geq 0.\n",
    "$$\n",
    "\n",
    "Therefore, $A^* A = R^* R$ always exists.\n",
    "\n",
    "Then the matrix $A R^{-1}$ is unitary:  \n",
    "\n",
    "$$\n",
    "   (A R^{-1})^* (AR^{-1})= R^{-*} A^* A R^{-1} = R^{-*} R^* R R^{-1} = I.\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Rank-deficient case\n",
    "When an $n \\times m$ matrix does not have full column rank, it is said to be **rank-deficient**.  \n",
    "\n",
    "The QR-decomposition, however, also exists.\n",
    "\n",
    "For any **rank-deficient matrix** there is a sequence of full-column rank matrices $A_k$ such that $A_k \\rightarrow A$.\n",
    "\n",
    "Each $A_k$ can be decomposed as $A_k = Q_k R_k$. \n",
    "\n",
    "The set of all unitary matrices is compact, thus there exists a converging subsequence $Q_{n_k} \\rightarrow Q$,\n",
    "\n",
    "and $Q^* A_k \\rightarrow Q^* A = R$, which is triangular."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR-decomposition and Gram matrices\n",
    "So, the simplest way to compute QR-decomposition is then\n",
    "\n",
    "$$A^* A = R^* R,$$\n",
    "\n",
    "(Cholesky factorization)\n",
    "\n",
    "$$Q = A R^{-1}.$$\n",
    "\n",
    "It is a **bad idea** for numerical stability. Let us do some demo (for a submatrix of the Hilbert matrix)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Built-in QR orth 8.77419464402e-16\n",
      "Via Gram matrix: 0.265317882364\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "n = 20\n",
    "r = 8\n",
    "#a = np.random.randn(n, r)\n",
    "a = [[1.0/(i+j+0.5) for i in xrange(r)] for j in xrange(n)]\n",
    "a = np.array(a)\n",
    "q, Rmat = np.linalg.qr(a)\n",
    "e = np.eye(r)\n",
    "print 'Built-in QR orth', np.linalg.norm(np.dot(q.T, q) - e)\n",
    "gram_matrix = a.T.dot(a)\n",
    "Rmat1 = np.linalg.cholesky(gram_matrix)\n",
    "q1 = np.dot(a, np.linalg.inv(Rmat1.T))\n",
    "print 'Via Gram matrix:', np.linalg.norm(np.dot(q1.T, q1) - e)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Second way: Gram-Schmidt orthogonalization\n",
    "QR-decomposition is a mathematical way of writing down the Gram-Schmidt orthogonalization process.  \n",
    "Given a sequence of vectors $a_1, \\ldots, a_m$ we want to find orthogonal basis $q_1, \\ldots, q_m$ such that every $a_i$ is a linear combination of such vectors.  \n",
    "\n",
    "**Gram-Schmidt:**\n",
    "1. $q_1 := a_1/\\Vert a_1 \\Vert$\n",
    "2. $q_2 := a_2 - (a_2, q_1) q_1, \\quad q_2 := q_2/\\Vert q_2 \\Vert$\n",
    "3. $q_3 := a_3 - (a_3, q_1) q_1 - (a_3, q_2) q_2, \\quad q_2 := q_3/\\Vert q_3 \\Vert$\n",
    "4. And go on  \n",
    "\n",
    "Note that the transformation from $Q$ to $A$ has triangular structure, since from the $k$-th vector we subtract only the previous ones."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Gram-Schmidt is unstable\n",
    "Gram-Schmidt can be very unstable (i.e., the produced vectors will be not orthogonal, especially if $q_k$ has small norm).  \n",
    "This is called **loss of orthogonality**.  \n",
    "\n",
    "There is a remedy, called **modified Gram-Schmidt** method. \n",
    "Instead of doing \n",
    "\n",
    "$$q_k := a_k - (a_k, q_1) q_1 - \\ldots - (a_k, q_{k-1}) q_{k-1}$$\n",
    "\n",
    "we do it step-by-step. First we set $q_k := a_k$ and orthogonalize sequentially:\n",
    "\n",
    "$$\n",
    "   q_k := q_k - (q_k, q_1)q_1, \\quad q_k := q_{k} - (q_k,q_2)q_2, \\ldots\n",
    "$$\n",
    "\n",
    "In exact arithmetic, it is the same. In floating point it is absolutely different!\n",
    "\n",
    "Note that the complexity is $\\mathcal{O}(nm^2)$ operations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR-decomposition: the (almost) practical way\n",
    "\n",
    "If $A = QR$, then  \n",
    "$$\n",
    "R = Q^* A,\n",
    "$$\n",
    "and we need to find a certain orthogonal matrix $Q$ that brings a matrix into orthogonal form.  \n",
    "For simplicity, we will look for an $n \\times n$ matrix such that\n",
    "$$\n",
    "   Q^* A = \\begin{bmatrix}\n",
    "   * & * & *  \\\\\n",
    "   0 & * & * \\\\\n",
    "   0 & 0 & * \\\\\n",
    "   &0_{(n-m) \\times m}\n",
    "   \\end{bmatrix}\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We will do it column-by-column.  \n",
    "First, we find a Householder matrix $H_1 = (I - 2 uu^{\\top})$ such that (we illustrate on a $4 \\times 3$ matrix)\n",
    "\n",
    "$$\n",
    "   H_1 A = \\begin{bmatrix}\n",
    "    * & * & * \\\\\n",
    "    0 & * & * \\\\\n",
    "    0 & * & * \\\\\n",
    "    0 & * & *\n",
    "   \\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Then, \n",
    "$$\n",
    "   H_2 H_1 A = \\begin{bmatrix}\n",
    "    * & * & * \\\\\n",
    "    0 & * & * \\\\\n",
    "    0 & 0 & * \\\\\n",
    "    0 & 0 & *\n",
    "   \\end{bmatrix},\n",
    "$$\n",
    "where\n",
    "$$\n",
    "  H_2 = \\begin{bmatrix}\n",
    "  1 & 0 \\\\\n",
    "  0 & H'_2, \n",
    "  \\end{bmatrix}\n",
    "$$\n",
    "and $H'_2$ is a $3 \\times 3$ Householder matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Finally, \n",
    "$$\n",
    "   H_3 H_2 H_1 A = \\begin{bmatrix}\n",
    "    * & * & * \\\\\n",
    "    0 & * & * \\\\\n",
    "    0 & 0 & * \\\\\n",
    "    0 & 0 & 0\n",
    "   \\end{bmatrix},\n",
    "$$\n",
    "\n",
    "You can try to implement by yourself, it is simple."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR-decomposition: real life\n",
    "\n",
    "In reality, since this is a dense matrix factorization, you should implement the algorithm in terms of blocks.  \n",
    "\n",
    "Instead of using Householder transformation, we use **block Householder** transformation of the form\n",
    "\n",
    "$$H = (I - 2UU^*), $$\n",
    "\n",
    "where $U^* U = I$.\n",
    "\n",
    "This allows us to use BLAS-3 operations. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Rank-revealing QR-decomposition\n",
    "\n",
    "The QR-decomposition can be also used to compute the (numerical) rank of the matrix.  \n",
    "\n",
    "It is done via so-called **rank-revealing** factorization.  \n",
    "\n",
    "It is based on the representation\n",
    "\n",
    "$$P A = Q R, $$\n",
    "\n",
    "where $P$ is the permutation matrix (it permutes columns), and  $R$ has the block form\n",
    "\n",
    "\n",
    "$$R = \\begin{bmatrix} R_{11} & R_{12} \\\\ 0 & R_{22}\\end{bmatrix},$$\n",
    "\n",
    "and the norm of $R_{22}$ is **small**, so you can find the numerical rank by looking at it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Summary\n",
    "\n",
    "LU and QR decompositions can be computed using **direct methods** in finite amount of operations.\n",
    "\n",
    "What about Schur form and SVD?  \n",
    "\n",
    "They can not be computed by **direct methods** they can only be computed by **iterative methods**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Schur form\n",
    "\n",
    "Every matrix can be written as a product\n",
    "\n",
    "$$A = U T U^*,$$\n",
    "\n",
    "with triangular $T$ and unitary $U$\n",
    "\n",
    "and this decomposition gives eigenvalues of the matrix.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR-algorithm\n",
    "The QR algorithm (Kublanovskaya and Francis independently proposed it in 1961).   \n",
    "\n",
    "\n",
    "\n",
    "<font color='red'> Do not **mix** QR algorithm and QR decomposition. </font>\n",
    "\n",
    "QR-decomposition is the representation of a matrix, whereas QR algorithm uses QR decomposition computes the eigenvalues!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Way to QR-algorithm\n",
    "\n",
    "Consider the equation\n",
    "\n",
    "$$A = Q T Q^*,$$\n",
    "\n",
    "and rewrite it in the form\n",
    "\n",
    "$$\n",
    "   Q T = A Q.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Derivation of the QR-algorithm as fixed-point\n",
    "Then we can write down the iterative process \n",
    "\n",
    "$$\n",
    "    Q_{k+1} R_{k+1} = A Q_k, \\quad Q_{k+1}^* A = R_{k+1} Q^*_k\n",
    "$$\n",
    "\n",
    "Introduce\n",
    "\n",
    "$$A_k = Q^* _k A Q_k = Q^*_k Q_{k+1} R_{k+1} = \\widehat{Q}_k R_{k+1}$$.\n",
    "\n",
    "and the new approximation reads\n",
    "\n",
    "$$A_{k+1} = Q^*_{k+1} A Q_{k+1} = \\widehat{Q}^*_k A_k \\widehat{Q}_k = R_{k+1} \\widehat{Q}_k.$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Final formulas\n",
    "The final formulas are then written in the form\n",
    "$$\n",
    "   A_0 = A, \\quad A_k = Q_k R_k, \\quad A_{k+1} = R_k Q_k.\n",
    "$$\n",
    "\n",
    "The matrices $A_{k}$ are unitary similar, and $A_k \\rightarrow T$.  \n",
    "\n",
    "The complexity of each step is $\\mathcal{O}(n^3)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "current a:\n",
      "[[  2.40047183e+00   1.43485636e-01   4.99605047e-03  -5.56291523e-05]\n",
      " [  1.43485636e-01   3.59286592e-01   1.60534687e-02  -1.94225770e-04]\n",
      " [  4.99605047e-03   1.60534687e-02   1.60692682e-02  -2.80885670e-04]\n",
      " [ -5.56291523e-05  -1.94225770e-04  -2.80885670e-04   2.40684296e-04]]\n",
      "current a:\n",
      "[[  2.41031150e+00   2.09437993e-02   3.18722176e-05   5.45279213e-09]\n",
      " [  2.09437993e-02   3.50196113e-01   6.87806955e-04   1.28079871e-07]\n",
      " [  3.18722176e-05   6.87806955e-04   1.53250849e-02   4.17732654e-06]\n",
      " [  5.45279205e-09   1.28079871e-07   4.17732654e-06   2.35678648e-04]]\n",
      "current a:\n",
      "[[  2.41051991e+00   3.04114601e-03   2.02621961e-07  -5.33227597e-13]\n",
      " [  3.04114601e-03   3.49989111e-01   3.00995528e-05  -8.62074652e-11]\n",
      " [  2.02621961e-07   3.00995528e-05   1.53236760e-02  -6.42429491e-08]\n",
      " [ -5.33147641e-13  -8.62074827e-11  -6.42429491e-08   2.35677492e-04]]\n",
      "current a:\n",
      "[[  2.41052431e+00   4.41545673e-04   1.28806658e-09   1.32059800e-16]\n",
      " [  4.41545673e-04   3.49984720e-01   1.31786000e-06   5.80333420e-14]\n",
      " [  1.28806664e-09   1.31786000e-06   1.53236733e-02   9.88053919e-10]\n",
      " [  5.21260158e-17   5.80510077e-14   9.88053906e-10   2.35677492e-04]]\n",
      "current a:\n",
      "[[  2.41052440e+00   6.41081268e-05   8.18816607e-12  -7.99356474e-17]\n",
      " [  6.41081268e-05   3.49984627e-01   5.77009674e-08  -2.14109052e-17]\n",
      " [  8.18822358e-12   5.77009674e-08   1.53236733e-02  -1.51962427e-11]\n",
      " [ -5.09637195e-21  -3.90911829e-17  -1.51962302e-11   2.35677492e-04]]\n",
      "current a:\n",
      "[[  2.41052440e+00   9.30787458e-06   5.19949264e-14   7.99300814e-17]\n",
      " [  9.30787458e-06   3.49984626e-01   2.52637036e-09  -1.76560777e-17]\n",
      " [  5.20524342e-14   2.52637031e-09   1.53236733e-02   2.33729939e-13]\n",
      " [  4.98273388e-25   2.63237618e-20   2.33717423e-13   2.35677492e-04]]\n",
      "current a:\n",
      "[[  2.41052440e+00   1.35141258e-06   2.73389069e-16  -7.99300126e-17]\n",
      " [  1.35141258e-06   3.49984625e-01   1.10614263e-10   1.76826923e-17]\n",
      " [  3.30896669e-16   1.10614211e-10   1.53236733e-02  -3.60708065e-15]\n",
      " [ -4.87162969e-29  -1.77262591e-23  -3.59456477e-15   2.35677492e-04]]\n",
      "current a:\n",
      "[[  2.41052440e+00   1.96211922e-07  -5.54040645e-17   7.99300027e-17]\n",
      " [  1.96211922e-07   3.49984625e-01   4.84316731e-12  -1.76827548e-17]\n",
      " [  2.10350596e-18   4.84311567e-12   1.53236733e-02   6.78001457e-17]\n",
      " [  4.76300288e-33   1.19367537e-26   5.52842647e-17   2.35677492e-04]]\n",
      "current a:\n",
      "[[  2.41052440e+00   2.84880569e-08  -5.74941943e-17  -7.99300013e-17]\n",
      " [  2.84880568e-08   3.49984625e-01   2.12101875e-13   1.76827614e-17]\n",
      " [  1.33719609e-20   2.12050235e-13   1.53236733e-02  -1.33661508e-17]\n",
      " [ -4.65679822e-37  -8.03813647e-30  -8.50269816e-19   2.35677492e-04]]\n",
      "current a:\n",
      "[[  2.41052440e+00   4.13618795e-09  -5.75074807e-17   7.99300011e-17]\n",
      " [  4.13618792e-09   3.49984625e-01   9.33601458e-15  -1.76827623e-17]\n",
      " [  8.50053870e-23   9.28437503e-15   1.53236733e-02   1.25289581e-17]\n",
      " [  4.55296169e-41   5.41283161e-33   1.30771163e-20   2.35677492e-04]]\n",
      "current a:\n",
      "[[  2.41052440e+00   6.00534165e-10  -5.75075651e-17  -7.99300010e-17]\n",
      " [  6.00534132e-10   3.49984625e-01   4.58145201e-16   1.76827624e-17]\n",
      " [  5.40378175e-25   4.06505655e-16   1.53236733e-02  -1.25160822e-17]\n",
      " [ -4.45144049e-45  -3.64496748e-36  -2.01125535e-22   2.35677492e-04]]\n",
      "current a:\n",
      "[[  2.41052440e+00   8.71917294e-11  -5.75075656e-17   7.99300010e-17]\n",
      " [  8.71916971e-11   3.49984625e-01   6.94379277e-17  -1.76827625e-17]\n",
      " [  3.43517725e-27   1.77983814e-17   1.53236733e-02   1.25158841e-17]\n",
      " [  4.35218299e-49   2.45449866e-39   3.09330281e-24   2.35677492e-04]]\n",
      "current a:\n",
      "[[  2.41052440e+00   1.26594160e-11  -5.75075656e-17  -7.99300010e-17]\n",
      " [  1.26593838e-11   3.49984625e-01   5.24188280e-17   1.76827625e-17]\n",
      " [  2.18373786e-29   7.79281605e-19   1.53236733e-02  -1.25158811e-17]\n",
      " [ -4.25513873e-53  -1.65284428e-42  -4.75748755e-26   2.35677492e-04]]\n",
      "current a:\n",
      "[[  2.41052440e+00   1.83805126e-12  -5.75075656e-17   7.99300010e-17]\n",
      " [  1.83801902e-12   3.49984625e-01   5.16736663e-17  -1.76827625e-17]\n",
      " [  1.38819941e-31   3.41199465e-20   1.53236733e-02   1.25158810e-17]\n",
      " [  4.16025833e-57   1.11301516e-45   7.31699713e-28   2.35677492e-04]]\n",
      "current a:\n",
      "[[  2.41052440e+00   2.66894672e-13  -5.75075656e-17  -7.99300010e-17]\n",
      " [  2.66862429e-13   3.49984625e-01   5.16410403e-17   1.76827625e-17]\n",
      " [  8.82476614e-34   1.49390251e-21   1.53236733e-02  -1.25158810e-17]\n",
      " [ -4.06749357e-61  -7.49497557e-49  -1.12535128e-29   2.35677492e-04]]\n",
      "current a:\n",
      "[[  2.41052440e+00   3.87780641e-14  -5.75075656e-17   7.99300010e-17]\n",
      " [  3.87458212e-14   3.49984625e-01   5.16396118e-17  -1.76827625e-17]\n",
      " [  5.60989272e-36   6.54087989e-23   1.53236733e-02   1.25158810e-17]\n",
      " [  3.97679726e-65   5.04707040e-52   1.73078584e-31   2.35677492e-04]]\n",
      "current a:\n",
      "[[  2.41052440e+00   5.65775813e-15  -5.75075656e-17  -7.99300010e-17]\n",
      " [  5.62551523e-15   3.49984625e-01   5.16395492e-17   1.76827625e-17]\n",
      " [  3.56620174e-38   2.86384884e-24   1.53236733e-02  -1.25158810e-17]\n",
      " [ -3.88812328e-69  -3.39866613e-55  -2.66194182e-33   2.35677492e-04]]\n",
      "current a:\n",
      "[[  2.41052440e+00   8.49012831e-16  -5.75075656e-17   7.99300010e-17]\n",
      " [  8.16769928e-16   3.49984625e-01   5.16395465e-17  -1.76827625e-17]\n",
      " [  2.26702996e-40   1.25390319e-25   1.53236733e-02   1.25158810e-17]\n",
      " [  3.80142654e-73   2.28864085e-58   4.09405604e-35   2.35677492e-04]]\n",
      "current a:\n",
      "[[  2.41052440e+00   1.50829928e-16  -5.75075656e-17  -7.99300010e-17]\n",
      " [  1.18587025e-16   3.49984625e-01   5.16395464e-17   1.76827625e-17]\n",
      " [  1.44114809e-42   5.49007053e-27   1.53236733e-02  -1.25158810e-17]\n",
      " [ -3.71666295e-77  -1.54115666e-61  -6.29664209e-37   2.35677492e-04]]\n",
      "current a:\n",
      "[[  2.41052440e+00   4.94605823e-17  -5.75075656e-17   7.99300010e-17]\n",
      " [  1.72176791e-17   3.49984625e-01   5.16395464e-17  -1.76827625e-17]\n",
      " [  9.16136026e-45   2.40376408e-28   1.53236733e-02   1.25158810e-17]\n",
      " [  3.63378941e-81   1.03780540e-64   9.68421076e-39   2.35677492e-04]]\n",
      "current a:\n",
      "[[  2.41052440e+00   3.47427422e-17  -5.75075656e-17  -7.99300010e-17]\n",
      " [  2.49983902e-18   3.49984625e-01   5.16395464e-17   1.76827625e-17]\n",
      " [  5.82386518e-47   1.05246038e-29   1.53236733e-02  -1.25158810e-17]\n",
      " [ -3.55276376e-85  -6.98851769e-68  -1.48942780e-40   2.35677492e-04]]\n",
      "current a:\n",
      "[[  2.41052440e+00   3.26058554e-17  -5.75075656e-17   7.99300010e-17]\n",
      " [  3.62952238e-19   3.49984625e-01   5.16395464e-17  -1.76827625e-17]\n",
      " [  3.70222376e-49   4.60807643e-31   1.53236733e-02   1.25158810e-17]\n",
      " [  3.47354481e-89   4.70602478e-71   2.29073410e-42   2.35677492e-04]]\n",
      "current a:\n",
      "[[  2.41052440e+00   3.22956003e-17  -5.75075656e-17  -7.99300010e-17]\n",
      " [  5.26971239e-20   3.49984625e-01   5.16395464e-17   1.76827625e-17]\n",
      " [  2.35349899e-51   2.01759313e-32   1.53236733e-02  -1.25158810e-17]\n",
      " [ -3.39609227e-93  -3.16900811e-74  -3.52314004e-44   2.35677492e-04]]\n",
      "current a:\n",
      "[[  2.41052440e+00   3.22505543e-17  -5.75075656e-17   7.99300010e-17]\n",
      " [  7.65110827e-21   3.49984625e-01   5.16395464e-17  -1.76827625e-17]\n",
      " [  1.49611634e-53   8.83379887e-34   1.53236733e-02   1.25158810e-17]\n",
      " [  3.32036676e-97   2.13399055e-77   5.41857552e-46   2.35677492e-04]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22440141e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  1.11086628e-021   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  9.51079275e-056   3.86777697e-035   1.53236733e-002  -1.25158810e-017]\n",
      " [ -3.24632976e-101  -1.43701610e-080  -8.33374783e-048   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22430645e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  1.61286945e-022   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  6.04599897e-058   1.69346155e-036   1.53236733e-002   1.25158810e-017]\n",
      " [  3.17394363e-105   9.67677792e-084   1.28172714e-049   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429266e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  2.34172909e-023   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  3.84343394e-060   7.41462611e-038   1.53236733e-002  -1.25158810e-017]\n",
      " [ -3.10317155e-109  -6.51628267e-087  -1.97129130e-051   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429066e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  3.39996218e-024   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  2.44326611e-062   3.24640855e-039   1.53236733e-002   1.25158810e-017]\n",
      " [  3.03397754e-113   4.38802463e-090   3.03183826e-053   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429037e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  4.93641339e-025   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  1.55318119e-064   1.42140255e-040   1.53236733e-002  -1.25158810e-017]\n",
      " [ -2.96632640e-117  -2.95486877e-093  -4.66295531e-055   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429033e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  7.16719063e-026   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  9.87355328e-067   6.22344715e-042   1.53236733e-002   1.25158810e-017]\n",
      " [  2.90018374e-121   1.98979043e-096   7.17160693e-057   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  1.04060616e-026   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  6.27660539e-069   2.72486457e-043   1.53236733e-002  -1.25158810e-017]\n",
      " [ -2.83551591e-125  -1.33991263e-099  -1.10299032e-058   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  1.51085862e-027   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  3.99003014e-071   1.19305053e-044   1.53236733e-002   1.25158810e-017]\n",
      " [  2.77229004e-129   9.02288913e-103   1.69639478e-060   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  2.19361931e-028   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  2.53645714e-073   5.22363419e-046   1.53236733e-002  -1.25158810e-017]\n",
      " [ -2.71047397e-133  -6.07595798e-106  -2.60904850e-062   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  3.18492122e-029   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  1.61242261e-075   2.28710800e-047   1.53236733e-002   1.25158810e-017]\n",
      " [  2.65003626e-137   4.09151270e-109   4.01270632e-064   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  4.62419489e-030   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  1.02501502e-077   1.00138386e-048   1.53236733e-002  -1.25158810e-017]\n",
      " [ -2.59094618e-141  -2.75519946e-112  -6.17152654e-066   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  6.71387984e-031   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  6.51600760e-080   4.38444378e-050   1.53236733e-002   1.25158810e-017]\n",
      " [  2.53317368e-145   1.85533435e-115   9.49178354e-068   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  9.74789851e-032   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  4.14221782e-082   1.91967815e-051   1.53236733e-002  -1.25158810e-017]\n",
      " [ -2.47668939e-149  -1.24937073e-118  -1.45983257e-069   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  1.41529976e-032   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  2.63320266e-084   8.40508943e-053   1.53236733e-002   1.25158810e-017]\n",
      " [  2.42146457e-153   8.41318556e-122   2.24521676e-071   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  2.05487718e-033   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  1.67392361e-086   3.68007149e-054   1.53236733e-002  -1.25158810e-017]\n",
      " [ -2.36747115e-157  -5.66538735e-125  -3.45313454e-073   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  2.98348118e-034   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  1.06411113e-088   1.61127687e-055   1.53236733e-002   1.25158810e-017]\n",
      " [  2.31468166e-161   3.81503696e-128   5.31090733e-075   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  4.33172361e-035   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  6.76454108e-091   7.05478996e-057   1.53236733e-002  -1.25158810e-017]\n",
      " [ -2.26306927e-165  -2.56902240e-131  -8.16815458e-077   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  6.28924007e-036   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  4.30021024e-093   3.08885843e-058   1.53236733e-002   1.25158810e-017]\n",
      " [  2.21260772e-169   1.72996387e-134   1.25625896e-078   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  9.13136299e-037   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  2.73363824e-095   1.35242105e-059   1.53236733e-002  -1.25158810e-017]\n",
      " [ -2.16327135e-173  -1.16494702e-137  -1.93212134e-080   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  1.32578482e-037   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  1.73777039e-097   5.92141961e-061   1.53236733e-002   1.25158810e-017]\n",
      " [  2.11503508e-177   7.84468147e-141   2.97159502e-082   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  1.92491021e-038   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  1.10469845e-099   2.59262530e-062   1.53236733e-002  -1.25158810e-017]\n",
      " [ -2.06787437e-181  -5.28256020e-144  -4.57030145e-084   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  2.79478183e-039   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  7.02255412e-102   1.13515109e-063   1.53236733e-002   1.25158810e-017]\n",
      " [  2.02176524e-185   3.55724352e-147   7.02910565e-086   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  4.05775056e-040   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  4.46422881e-104   4.97012816e-065   1.53236733e-002  -1.25158810e-017]\n",
      " [ -1.97668425e-189  -2.39542588e-150  -1.08107368e-087   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  5.89145793e-041   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  2.83790463e-106   2.17611330e-066   1.53236733e-002   1.25158810e-017]\n",
      " [  1.93260847e-193   1.61306504e-153   1.66268707e-089   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  8.55382214e-042   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  1.80405240e-108   9.52786117e-068   1.53236733e-002  -1.25158810e-017]\n",
      " [ -1.88951548e-197  -1.08622807e-156  -2.55720618e-091   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  1.24193152e-042   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  1.14683384e-110   4.17166415e-069   1.53236733e-002   1.25158810e-017]\n",
      " [  1.84738337e-201   7.31459291e-160   3.93297304e-093   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  1.80316341e-043   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  7.29040829e-113   1.82651504e-070   1.53236733e-002  -1.25158810e-017]\n",
      " [ -1.80619071e-205  -4.92560183e-163  -6.04889706e-095   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  2.61801735e-044   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  4.63450337e-115   7.99718550e-072   1.53236733e-002   1.25158810e-017]\n",
      " [  1.76591657e-209   3.31686995e-166   9.30317988e-097   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  3.80110577e-045   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  2.94614796e-117   3.50147546e-073   1.53236733e-002  -1.25158810e-017]\n",
      " [ -1.72654045e-213  -2.23355980e-169  -1.43082541e-098   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  5.51883474e-046   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  1.87286255e-119   1.53308065e-074   1.53236733e-002   1.25158810e-017]\n",
      " [  1.68804233e-217   1.50406542e-172   2.20060385e-100   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  8.01280961e-047   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  1.19057636e-121   6.71241686e-076   1.53236733e-002  -1.25158810e-017]\n",
      " [ -1.65040264e-221  -1.01282839e-175  -3.38452006e-102   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  1.16338178e-047   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  7.56847895e-124   2.93895432e-077   1.53236733e-002   1.25158810e-017]\n",
      " [  1.61360223e-225   6.82032402e-179   5.20537853e-104   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  1.68911685e-048   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  4.81127254e-126   1.28678726e-078   1.53236733e-002  -1.25158810e-017]\n",
      " [ -1.57762239e-229  -4.59276420e-182  -8.00585170e-106   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  2.45243287e-049   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  3.05851989e-128   5.63404964e-080   1.53236733e-002   1.25158810e-017]\n",
      " [  1.54244482e-233   3.09273913e-185   1.23129684e-107   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  3.56069326e-050   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  1.94429725e-130   2.46680367e-081   1.53236733e-002  -1.25158810e-017]\n",
      " [ -1.50805164e-237  -2.08263149e-188  -1.89372970e-109   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  5.16977923e-051   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  1.23598731e-132   1.08006154e-082   1.53236733e-002   1.25158810e-017]\n",
      " [  1.47442534e-241   1.40243122e-191   2.91254883e-111   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  7.50601507e-052   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  7.85715578e-135   4.72892492e-084   1.53236733e-002  -1.25158810e-017]\n",
      " [ -1.44154885e-245  -9.44388550e-195  -4.47948864e-113   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  1.08980016e-052   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  4.99478403e-137   2.07050525e-085   1.53236733e-002   1.25158810e-017]\n",
      " [  1.40940543e-249   6.35945435e-198   6.88943590e-115   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  1.58228351e-053   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  3.17517792e-139   9.06546848e-087   1.53236733e-002  -1.25158810e-017]\n",
      " [ -1.37797873e-253  -4.28241740e-201  -1.05959253e-116   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  2.29732129e-054   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  2.01845661e-141   3.96921084e-088   1.53236733e-002   1.25158810e-017]\n",
      " [  1.34725279e-257   2.88375351e-204   1.62964915e-118   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  3.33548638e-055   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  1.28313032e-143   1.73787320e-089   1.53236733e-002  -1.25158810e-017]\n",
      " [ -1.31721196e-261  -1.94190186e-207  -2.50639397e-120   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  4.84280080e-056   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  8.15684331e-146   7.60907740e-091   1.53236733e-002   1.25158810e-017]\n",
      " [  1.28784099e-265   1.30766476e-210   3.85482407e-122   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  7.03127430e-057   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  5.18529502e-148   3.33154680e-092   1.53236733e-002  -1.25158810e-017]\n",
      " [ -1.25912492e-269  -8.80573395e-214  -5.92870426e-124   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  1.02087243e-057   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  3.29628552e-150   1.45867935e-093   1.53236733e-002   1.25158810e-017]\n",
      " [  1.23104915e-273   5.92972702e-217   9.11832382e-126   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  1.48220717e-058   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  2.09544455e-152   6.38665933e-095   1.53236733e-002  -1.25158810e-017]\n",
      " [ -1.20359942e-277  -3.99304167e-220  -1.40239462e-127   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  2.15202021e-059   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  1.33207146e-154   2.79632515e-096   1.53236733e-002   1.25158810e-017]\n",
      " [  1.17676176e-281   2.68888967e-223   2.15687740e-129   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  3.12452339e-060   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  8.46796153e-157   1.22433872e-097   1.53236733e-002  -1.25158810e-017]\n",
      " [ -1.15052251e-285  -1.81068175e-226  -3.31726896e-131   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  4.53650312e-061   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  5.38307248e-159   5.36062591e-099   1.53236733e-002   1.25158810e-017]\n",
      " [  1.12486835e-289   1.21930194e-229   5.10194661e-133   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  6.58655995e-062   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  3.42201240e-161   2.34708824e-100   1.53236733e-002  -1.25158810e-017]\n",
      " [ -1.09978622e-293  -8.21070420e-233  -7.84677382e-135   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  9.56304245e-063   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  2.17536897e-163   1.02764552e-101   1.53236733e-002   1.25158810e-017]\n",
      " [  1.07526337e-297   5.52903765e-236   1.20683073e-136   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  1.38846047e-063   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  1.38287932e-165   4.49942742e-103   1.53236733e-002  -1.25158810e-017]\n",
      " [ -1.05128732e-301  -3.72321991e-239  -1.85610091e-138   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  2.01590914e-064   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  8.79094641e-168   1.97002241e-104   1.53236733e-002   1.25158810e-017]\n",
      " [  1.02784588e-305   2.50719337e-242   2.85467590e-140   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  2.92690339e-065   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  5.58839357e-170   8.62551598e-106   1.53236733e-002  -1.25158810e-017]\n",
      " [ -1.00492714e-309  -1.68832857e-245  -4.39048030e-142   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  4.24957817e-066   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  3.55253475e-172   3.77658271e-107   1.53236733e-002   1.25158810e-017]\n",
      " [  9.82519441e-314   1.13691007e-248   6.75254142e-144   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  6.16997291e-067   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  2.25834187e-174   1.65353319e-108   1.53236733e-002  -1.25158810e-017]\n",
      " [ -9.60611341e-318  -7.65588239e-252  -1.03853821e-145   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  8.95819872e-068   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  1.43562508e-176   7.23980441e-110   1.53236733e-002   1.25158810e-017]\n",
      " [  9.38724727e-322   5.15542406e-255   1.59726769e-147   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  1.30064306e-068   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  9.12625058e-179   3.16986488e-111   1.53236733e-002  -1.25158810e-017]\n",
      " [  0.00000000e+000  -3.47163082e-258  -2.45659142e-149   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  1.88840683e-069   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  5.80154600e-181   1.38788879e-112   1.53236733e-002   1.25158810e-017]\n",
      " [  0.00000000e+000   2.33777482e-261   3.77822794e-151   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  2.74178247e-070   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  3.68803549e-183   6.07671101e-114   1.53236733e-002  -1.25158810e-017]\n",
      " [  0.00000000e+000  -1.57424317e-264  -5.81089971e-153   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  3.98080066e-071   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  2.34447952e-185   2.66061785e-115   1.53236733e-002   1.25158810e-017]\n",
      " [  0.00000000e+000   1.06008566e-267   8.93714089e-155   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  5.77973419e-072   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  1.49038268e-187   1.16492084e-116   1.53236733e-002  -1.25158810e-017]\n",
      " [  0.00000000e+000  -7.13855155e-271  -1.37452875e-156   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  8.39161016e-073   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  9.47434391e-190   5.10047159e-118   1.53236733e-002   1.25158810e-017]\n",
      " [  0.00000000e+000   4.80705666e-274   2.11401981e-158   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  1.21837993e-073   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  6.02282850e-192   2.23318267e-119   1.53236733e-002  -1.25158810e-017]\n",
      " [  0.00000000e+000  -3.23704236e-277  -3.25135416e-160   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  1.76896879e-074   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  3.82870449e-194   9.77773282e-121   1.53236733e-002   1.25158810e-017]\n",
      " [  0.00000000e+000   2.17980439e-280   5.00056990e-162   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  2.56837009e-075   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  2.43390263e-196   4.28106757e-122   1.53236733e-002  -1.25158810e-017]\n",
      " [  0.00000000e+000  -1.46786686e-283  -7.69085684e-164   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  3.72902280e-076   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  1.54722883e-198   1.87441607e-123   1.53236733e-002   1.25158810e-017]\n",
      " [  0.00000000e+000   9.88452503e-287   1.18285076e-165   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  5.41417730e-077   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  9.83571422e-201   8.20691462e-125   1.53236733e-002  -1.25158810e-017]\n",
      " [  0.00000000e+000  -6.65617829e-290  -1.81921981e-167   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  7.86085723e-078   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  6.25255115e-203   3.59330293e-126   1.53236733e-002   1.25158810e-017]\n",
      " [  0.00000000e+000   4.48222948e-293   2.79795291e-169   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  1.14131978e-078   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  3.97473889e-205   1.57328626e-127   1.53236733e-002  -1.25158810e-017]\n",
      " [  0.00000000e+000  -3.01830573e-296  -4.30324057e-171   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  1.65708497e-079   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  2.52673651e-207   6.88845246e-129   1.53236733e-002   1.25158810e-017]\n",
      " [  0.00000000e+000   2.03250850e-299   6.61836707e-173   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  2.40592571e-080   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  1.60624322e-209   3.01602948e-130   1.53236733e-002  -1.25158810e-017]\n",
      " [  0.00000000e+000  -1.36867871e-302  -1.01790225e-174   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  3.49316940e-081   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  1.02108679e-211   1.32053373e-131   1.53236733e-002   1.25158810e-017]\n",
      " [  0.00000000e+000   9.21659814e-306   1.56552966e-176   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  5.07174117e-082   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  6.49103587e-214   5.78180466e-133   1.53236733e-002  -1.25158810e-017]\n",
      " [  0.00000000e+000  -6.20640044e-309  -2.40777845e-178   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  7.36367338e-083   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  4.12634333e-216   2.53149650e-134   1.53236733e-002   1.25158810e-017]\n",
      " [  0.00000000e+000   4.17935184e-312   3.70315379e-180   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017  -7.99300010e-017]\n",
      " [  1.06913353e-083   3.49984625e-001   5.16395464e-017   1.76827625e-017]\n",
      " [  2.62311126e-218   1.10838655e-135   1.53236733e-002  -1.25158810e-017]\n",
      " [  0.00000000e+000  -2.81434980e-315  -5.69543596e-182   2.35677492e-004]]\n",
      "current a:\n",
      "[[  2.41052440e+000   3.22429032e-017  -5.75075656e-017   7.99300010e-017]\n",
      " [  1.55227759e-084   3.49984625e-001   5.16395464e-017  -1.76827625e-017]\n",
      " [  1.66750853e-220   4.85294271e-137   1.53236733e-002   1.25158810e-017]\n",
      " [  0.00000000e+000   1.89516665e-318   8.75955810e-184   2.35677492e-004]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "n = 4\n",
    "a = [[1.0/(i + j + 0.5) for i in xrange(n)] for j in xrange(n)]\n",
    "niters = 100\n",
    "for k in xrange(niters):\n",
    "    q, rmat = np.linalg.qr(a)\n",
    "    a = rmat.dot(q)\n",
    "    print 'current a:'\n",
    "    print a\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Convergence and complexity of the QR-algorithm\n",
    "The convergence of the QR-algorithm is from **largest eigenvalues** to the smallest.\n",
    "\n",
    "At least 2-3 iterations is needed for an eigenvalue. \n",
    "\n",
    "Each step is one QR-factorization and one matrix-by-matrix product. $\\mathcal{O}(n^3)$ complexity.\n",
    "\n",
    "It means, $\\mathcal{O}(n^4)$ complexity? \n",
    "\n",
    "Fortunately, not. \n",
    "\n",
    "We can also speedup the QR-algorithm by using **shifts**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## A few words about the SVD\n",
    "\n",
    "And the last, but not the least, is the **singular value decomposition** of matrices.\n",
    "\n",
    "$$A = U \\Sigma V^*.$$\n",
    "\n",
    "We can compute via eigendecomposition of \n",
    "\n",
    "$$A^* A = U^* \\Sigma^2 U,$$\n",
    "\n",
    "but it is a **bad idea** (c.f. Gram matrices)\n",
    "\n",
    "We will discuss methods for computing SVD later (there will a special SVD lecture).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Summary of todays lecture\n",
    "\n",
    "- LU decomposition and Gaussian elimination\n",
    "- QR decomposition and Gram-Schmidt algorithm, reduction to a simpler form by Householder transformations\n",
    "- Schur decomposition and QR-algorithm\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Next lecture\n",
    "- No Test\n",
    "- SVD: how to compute it, what are the applications, computing best rank-$r$ approximation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##### Questions?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>\n",
       "<style>\n",
       "    @font-face {\n",
       "        font-family: \"Computer Modern\";\n",
       "        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "    }\n",
       "    div.cell{\n",
       "        /*width:80%;*/\n",
       "        /*margin-left:auto !important;\n",
       "        margin-right:auto;*/\n",
       "    }\n",
       "    h1 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\n",
       "    h2 {\n",
       "        font-family: 'Fenix', serif;\n",
       "    }\n",
       "    h3{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "        margin-top:12px;\n",
       "        margin-bottom: 3px;\n",
       "       }\n",
       "\th4{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "       }\n",
       "    h5 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\t   \n",
       "    div.text_cell_render{\n",
       "        font-family: 'Alegreya Sans',Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
       "        line-height: 1.2;\n",
       "        font-size: 120%;\n",
       "        /*width:70%;*/\n",
       "        /*margin-left:auto;*/\n",
       "        margin-right:auto;\n",
       "    }\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "    .text_cell_render h1 {\n",
       "        font-weight: 200;\n",
       "        font-size: 50pt;\n",
       "\t\tline-height: 110%;\n",
       "        color:#CD2305;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\t\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 300;\n",
       "        font-size: 16pt;\n",
       "        color: #CD2305;\n",
       "        font-style: italic;\n",
       "        margin-bottom: .5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\n",
       "    \n",
       "    li {\n",
       "        line-height: 110%;\n",
       "    }\n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"./styles/custom.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}