{
"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 positive definite 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",
" Do not **mix** QR algorithm and QR decomposition. \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": [
"\n",
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"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
}