section-matrix-factor.ipynb

{
"cells": [
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":["%%html\n<link href=\"https:\/\/pretextbook.org\/beta\/mathbook-content.css\" rel=\"stylesheet\" type=\"text\/css\" \/>\n<link href=\"https:\/\/aimath.org\/mathbook\/mathbook-add-on.css\" rel=\"stylesheet\" type=\"text\/css\" \/>\n<link href=\"https:\/\/fonts.googleapis.com\/css?family=Open+Sans:400,400italic,600,600italic\" rel=\"stylesheet\" type=\"text\/css\" \/>\n<link href=\"https:\/\/fonts.googleapis.com\/css?family=Inconsolata:400,700&subset=latin,latin-ext\" rel=\"stylesheet\" type=\"text\/css\" \/><!-- Hide this cell. -->\n<script>\nvar cell = $(\".container .cell\").eq(0), ia = cell.find(\".input_area\")\nif (cell.find(\".toggle-button\").length == 0) {\nia.after(\n    $('<button class=\"toggle-button\">Toggle hidden code<\/button>').click(\n        function (){ ia.toggle() }\n        )\n    )\nia.hide()\n}\n<\/script>\n"], "outputs":[]},
{"cell_type":"markdown", "metadata":{}, "source":["**Important:** to view this notebook properly you will need to execute the cell above, which assumes you have an Internet connection.  It should already be selected, or place your cursor anywhere above to select.  Then press the \"Run\" button in the menu bar above (the right-pointing arrowhead), or press Shift-Enter on your keyboard."]},
{"cell_type":"markdown", "metadata":{}, "source":["$\\newcommand{\\spn}{\\operatorname{span}}\n\\newcommand{\\bbm}{\\begin{bmatrix}}\n\\newcommand{\\ebm}{\\end{bmatrix}}\n\\newcommand{\\R}{\\mathbb{R}}\n\\renewcommand{\\C}{\\mathbb{C}}\n\\newcommand{\\im}{\\operatorname{im}}\n\\newcommand{\\nll}{\\operatorname{null}}\n\\newcommand{\\csp}{\\operatorname{col}}\n\\newcommand{\\rank}{\\operatorname{rank}}\n\\newcommand{\\diag}{\\operatorname{diag}}\n\\newcommand{\\tr}{\\operatorname{tr}}\n\\newcommand{\\dotp}{\\!\\boldsymbol{\\cdot}\\!}\n\\newcommand{\\len}[1]{\\lVert #1\\rVert}\n\\newcommand{\\abs}[1]{\\lvert #1\\rvert}\n\\newcommand{\\proj}[2]{\\operatorname{proj}_{#1}{#2}}\n\\newcommand{\\bz}{\\overline{z}}\n\\newcommand{\\zz}{\\mathbf{z}}\n\\newcommand{\\uu}{\\mathbf{u}}\n\\newcommand{\\vv}{\\mathbf{v}}\n\\newcommand{\\ww}{\\mathbf{w}}\n\\newcommand{\\xx}{\\mathbf{x}}\n\\newcommand{\\yy}{\\mathbf{y}}\n\\newcommand{\\zer}{\\mathbf{0}}\n\\newcommand{\\vecq}{\\mathbf{q}}\n\\newcommand{\\vecp}{\\mathbf{p}}\n\\newcommand{\\vece}{\\mathbf{e}}\n\\newcommand{\\basis}[2]{\\{\\mathbf{#1}_1,\\mathbf{#1}_2,\\ldots,\\mathbf{#1}_{#2}\\}}\n\\newcommand{\\gt}{>}\n\\newcommand{\\amp}{&}\n$"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><h6 class=\"heading hide-type\"><span xmlns:svg=\"http:\/\/www.w3.org\/2000\/svg\" class=\"type\">Section<\/span> <span class=\"codenumber\">4.5<\/span> <span class=\"title\">Matrix Factorizations and Eigenvalues<\/span><\/h6><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-691\">This section is a rather rapid tour of some cool ideas that get a lot of use in applied linear algebra. We are rather light on details here. The interested reader can consult sections 8.3–8.6 in Nicholson.<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><h6 class=\"heading hide-type\"><span xmlns:svg=\"http:\/\/www.w3.org\/2000\/svg\" class=\"type\">Subsection<\/span> <span class=\"codenumber\">4.5.1<\/span> <span class=\"title\">Matrix Factorizations<\/span><\/h6><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-692\">Recall that an $n\\times n$ matrix $A$ is symmetric if $A^T=A$ and hermitian if $A^H=A\\text{,}$ where $A^H$ is the conjugate transpose of a complex matrix. In either case, the corresponding matrix transformation $T_A$ is said to be <dfn xmlns:svg=\"http:\/\/www.w3.org\/2000\/svg\" class=\"terminology\">self-adjoint<\/dfn>, which means that it satisfies the condition<\/p><div class=\"displaymath\">\n\\begin{equation*}\n\\langle u,T_Av\\rangle = \\langle T_Au,v\\rangle\n\\end{equation*}\n<\/div><p class=\"continuation\">for all $u,v\\in \\mathbb{R}^n$ (or $\\mathbb{C}^n$).<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-693\">All such matrices (or operators) can be diagonalized, in the sense that there is an orthonormal basis of eigenvectors for that matrix. These eigenvectors can be arranged to form an orthogonal matrix $P$ (or unitary matrix $U$) such that<\/p><div xmlns:svg=\"http:\/\/www.w3.org\/2000\/svg\" class=\"displaymath\">\n\\begin{equation*}\nP^TAP = D \\quad \\text{ (or } U^HAU=D)\\text{,}\n\\end{equation*}\n<\/div><p class=\"continuation\">where $D$ is a diagonal matrix whose entries are the eigenvalues of $A\\text{.}$<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><h6 class=\"heading\"><span class=\"title\">Positive Operators.<\/span><\/h6><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><article class=\"definition definition-like\" id=\"def-positive-op\"><h6 xmlns:svg=\"http:\/\/www.w3.org\/2000\/svg\" class=\"heading\"><span class=\"type\">Definition<\/span><span class=\"space\"> <\/span><span class=\"codenumber\">4.5.1<\/span><span class=\"period\">.<\/span><\/h6><p id=\"p-694\">A symmetric\/hermitian matrix $A$ (or operator $T$) is <dfn class=\"terminology\">positive<\/dfn> if $\\xx^TA\\xx\\geq 0$ ($\\langle \\xx,T\\xx\\rangle\\geq 0$) for all vectors $\\xx\\neq \\zer\\text{.}$ It is <dfn class=\"terminology\">positive-definite<\/dfn> if $\\xx^TA\\xx\\gt 0$ for all nonzero $\\xx\\text{.}$<\/p><\/article><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><aside class=\"aside aside-like\" id=\"aside-3\"><p id=\"p-695\">Some books will define positive-definite operators by the condition $\\xx^TA\\xx$ without the requirement that $A$ must be symmetric\/hermitian. However, we will stick to the simpler definition.<\/p><\/aside><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-696\">This is equivalent to requiring that all the eigenvalues of $A$ are non-negative. Every positive matrix $A$ has a unique positive square root: a matrix $R$ such that $R^2=A\\text{.}$<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><article class=\"theorem theorem-like\" id=\"thm-positive-prod\"><h6 xmlns:svg=\"http:\/\/www.w3.org\/2000\/svg\" class=\"heading\"><span class=\"type\">Theorem<\/span><span class=\"space\"> <\/span><span class=\"codenumber\">4.5.2<\/span><span class=\"period\">.<\/span><\/h6><p id=\"p-697\">For any $n\\times n$ matrix $U\\text{,}$ the matrix $A=U^TU$ is positive. Moreover, if $U$ is invertible, then $A$ is positive-definite.<\/p><\/article><article class=\"hiddenproof\" id=\"proof-25\"><h6 class=\"heading\"><span class=\"type\">Proof<span class=\"period\">.<\/span><\/span><\/h6><p id=\"p-698\">For any $\\xx\\neq \\zer$ in $\\R^n\\text{,}$<\/p><div class=\"displaymath\">\n\\begin{equation*}\n\\xx^T A\\xx = \\xx^TU^T U\\xx = (U\\xx)^T(U\\xx) = \\len{U\\xx}^2\\geq 0\\text{.}\n\\end{equation*}\n<\/div><\/article><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-699\">What is interesting is that the converse to the above statement is also true. The <dfn xmlns:svg=\"http:\/\/www.w3.org\/2000\/svg\" class=\"terminology\">Cholesky factorization<\/dfn> of a positive-definite matrix $A$ is given by $A=U^TU\\text{,}$ where $U$ is upper-triangular, with positive diagonal entries. (See Nicholson for details.)<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-700\">Even better is that there is a very simple algorithm for obtaining the factorization: Carry the matrix $A$ to triangular form, using only row operations of the type $R_i+kR_j\\to R_i\\text{.}$ Then divide each row by the square root of the diagonal entry.<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-701\">The SymPy library contains the <code class=\"code-inline tex2jax_ignore\">cholesky()<\/code> algorithm. Note however that it produces a lower triangular matrix, rather than upper triangular. (That is, the output gives $L=U^T$ rather than $U\\text{,}$ so you will have $A=LL^T\\text{.}$) Let's give it a try. First, enter a positive-definite matrix. (We'll try the one from Example 8.3.3 in Nicholson.)<\/p><\/div>"]},
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":["from sympy import Matrix,init_printing\ninit_printing()\nA = Matrix([[10,5,2],[5,3,2],[2,2,3]])\nA"], "outputs":[]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-702\">Next, find the Cholesky factorization:<\/p><\/div>"]},
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":["L = A.cholesky()\nL, L*L.T"], "outputs":[]},
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":["L*L.T == A"], "outputs":[]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><h6 class=\"heading\"><span class=\"title\">Singular Value Decomposition.<\/span><\/h6><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-703\">For any $n\\times n$ matrix $A\\text{,}$ the matrices $A^TA$ and $AA^T$ are both positive. (Exercise!) This means that we can define $\\sqrt{A^TA}\\text{,}$ even if $A$ itself is not symmetric or positive.<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><ul id=\"p-704\" class=\"disc\"><li id=\"li-129\"><p id=\"p-705\">Since $A^TA$ is symmetric, we know that it can be diagonalized.<\/p><\/li><li id=\"li-130\"><p id=\"p-706\">Since $A^TA$ is positive, we know its eigenvalues are non-negative.<\/p><\/li><li id=\"li-131\"><p id=\"p-707\">This means we can define the <dfn xmlns:svg=\"http:\/\/www.w3.org\/2000\/svg\" class=\"terminology\">singular values<\/dfn> $\\sigma_i = \\sqrt{\\lambda_i}$ for each $i=1,\\ldots, n\\text{.}$<\/p><\/li><li id=\"li-132\"><p id=\"p-708\"><em class=\"alert\">Note:<\/em> it's possible to do this even if $A$ is not a square matrix!<\/p><\/li><\/ul><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-709\">The SymPy library has a function for computing the singular values of a matrix. Given a matrix <code class=\"code-inline tex2jax_ignore\">A<\/code>, the command <code class=\"code-inline tex2jax_ignore\">A.singular_values()<\/code> will return its singular values. Try this for a few different matrices below:<\/p><\/div>"]},
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":["A = Matrix([[1,2,3],[4,5,6]])\nA.singular_values()"], "outputs":[]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-710\">In fact, SymPy can even return singular values for a matrix with variable entries! Try the following example from the <a class=\"external\" href=\"https:\/\/docs.sympy.org\/latest\/modules\/matrices\/matrices.html#sympy.matrices.matrices.MatrixEigen.singular%5Fvalues\" target=\"_blank\">SymPy documentation<\/a>.<\/p><\/div>"]},
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":["from sympy import Symbol\nx = Symbol('x', real=True)\nM = Matrix([[0,1,0],[0,x,0],[-1,0,0]])\nM,M.singular_values()"], "outputs":[]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-711\">For an $n\\times n$ matrix $A\\text{,}$ we might not be able to diagonalize $A$ (with a single orthonormal basis). However, it turns out that it's <em xmlns:svg=\"http:\/\/www.w3.org\/2000\/svg\" class=\"emphasis\">always<\/em> possible to find a pair of orthonormal bases $\\{e_1,\\ldots, e_n\\}, \\{f_1,\\ldots, f_n\\}$ such that<\/p><div class=\"displaymath\">\n\\begin{equation*}\nAx = \\sigma_1(x\\cdot e_1)f_1+\\cdots + \\sigma_n(x\\cdot e_n)f_n\\text{.}\n\\end{equation*}\n<\/div><p class=\"continuation\">In matrix form, $A = P\\Sigma_A Q^T$ for orthogonal matrices $P,Q\\text{.}$<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-712\">In fact, this can be done even if $A$ is not square, which is arguably the more interesting case! Let $A$ be an $m\\times n$ matrix. We will find an $m\\times m$ orthogonal matrix $P$ and $n\\times n$ orthogonal matrix $Q\\text{,}$ such that $A=P\\Sigma_A Q^T\\text{,}$ where $\\Sigma_A$ is also $m\\times n\\text{.}$<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><aside class=\"aside aside-like\" id=\"aside-4\"><p id=\"p-713\">If $A$ is symmetric and positive-definite, the singular values of $A$ are just the eigenvalues of $A\\text{,}$ and the singular value decomposition is the same as diagonalization.<\/p><\/aside><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-714\">The basis $\\{f_1,\\ldots, f_n\\}$ is an orthonormal basis for $A^TA\\text{,}$ and the matrix $Q$ is the matrix whose columns are the vectors $f_i\\text{.}$ As a result, $Q$ is orthogonal.<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-715\">The matrix $\\Sigma_A$ is the same size as $A\\text{.}$ First, we list the positive singular values of $A$ in decreasing order:<\/p><div xmlns:svg=\"http:\/\/www.w3.org\/2000\/svg\" class=\"displaymath\">\n\\begin{equation*}\n\\sigma_1\\geq \\sigma_2\\geq \\cdots \\geq \\sigma_k\\gt 0\\text{.}\n\\end{equation*}\n<\/div><p class=\"continuation\">Then, we let $D_A = \\operatorname{diag}(\\sigma_1,\\ldots, \\sigma_k)\\text{,}$ and set<\/p><div class=\"displaymath\">\n\\begin{equation*}\n\\Sigma_A = \\begin{bmatrix}D_A\\amp 0\\\\0\\amp 0\\end{bmatrix}\\text{.}\n\\end{equation*}\n<\/div><p class=\"continuation\">That is, we put $D_A$ in the upper-left, and then fill in zeros as needed, until $\\Sigma_A$ is the same size as $A\\text{.}$<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-716\">Next, we compute the vectors $e_i = \\frac{1}{\\len{Af_i}}Af_i\\text{,}$ for $i=1,\\ldots, k\\text{.}$ As shown in Nicolson, $\\{e_1,\\ldots, e_r\\}$ will be an orthonormal basis for the column space of $A\\text{.}$ The matrix $P$ is constructed by extending this to an orthonormal basis of $\\R^m\\text{.}$<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-717\">All of this is a lot of work to do by hand, but it turns out that it can be done numerically, and more importantly, <em xmlns:svg=\"http:\/\/www.w3.org\/2000\/svg\" class=\"emphasis\">efficiently<\/em>, by a computer. The SymPy library does not have an <abbr class=\"initialism\">SVD<\/abbr> algorithm, but the <code class=\"code-inline tex2jax_ignore\">mpmath<\/code> library does, and it works well with SymPy.<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-718\">Example 8.6.1 in Nicolson shows that for the matrix $A = \\begin{bmatrix}1\\amp 0\\amp 1\\\\-1\\amp 1\\amp 0\\end{bmatrix}$ we should have<\/p><div xmlns:svg=\"http:\/\/www.w3.org\/2000\/svg\" class=\"displaymath\">\n\\begin{equation*}\nP = \\frac{1}{\\sqrt{2}}\\bbm 1\\amp 1\\\\-1\\amp 1\\ebm, \\Sigma_A = \\bbm \\sqrt{3}\\amp 0\\amp 0\\\\0\\amp 1\\amp 0\\ebm, Q = \\frac{1}{\\sqrt{6}}\\bbm2\\amp -1\\amp 1\\\\0\\amp \\sqrt{3}\\amp \\sqrt{3}\\\\-\\sqrt{2}\\amp -\\sqrt{2}\\amp \\sqrt{2}\\ebm\\text{.}\n\\end{equation*}\n<\/div><p class=\"continuation\">Let us test this on the computer. First, we import the <code class=\"code-inline tex2jax_ignore\">mpmath<\/code> library, and define $A\\text{.}$ (Importantly, the <code class=\"code-inline tex2jax_ignore\">svd<\/code> routine from <code class=\"code-inline tex2jax_ignore\">mpmath<\/code> will accept a SymPy matrix as input. We do not, for example, have to first convert it to a NumPy array.) Since both SymPy and mpmath have similar functions, we <code class=\"code-inline tex2jax_ignore\">import as<\/code>  for both libraries, so we can distinguish between them as needed.<\/p><\/div>"]},
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":["import sympy as sy\nimport mpmath as mp\nsy.init_printing()\nmp.pretty = True\nA = mp.matrix([[1,0,1],[-1,1,0]])\nA"], "outputs":[]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-719\">Note that we defined <code class=\"code-inline tex2jax_ignore\">A<\/code> as a mpmath <code class=\"code-inline tex2jax_ignore\">matrix<\/code>. The <code class=\"code-inline tex2jax_ignore\">svd<\/code> command from the mpmath library will also work on a SymPy <code class=\"code-inline tex2jax_ignore\">Matrix<\/code>, but certain mpmath commands, like the <code class=\"code-inline tex2jax_ignore\">chop<\/code> command used to round decimals, will not. Next, we apply the <code class=\"code-inline tex2jax_ignore\">svd<\/code> algorithm.<\/p><\/div>"]},
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":["P,S,QT = mp.svd(A)\nP,S,QT"], "outputs":[]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-720\">Note that the input isn't quite so nice this time: the algorithm is numerical. For some reason, if you define <code class=\"code-inline tex2jax_ignore\">A<\/code> as a SymPy matrix, then <code class=\"code-inline tex2jax_ignore\">P<\/code> will be a SymPy matrix, but <code class=\"code-inline tex2jax_ignore\">S<\/code> and <code class=\"code-inline tex2jax_ignore\">Q<\/code> will not. (I have no idea why.) Also, the matrix <code class=\"code-inline tex2jax_ignore\">S<\/code> is a column matrix, containing the positive singular values of $A\\text{.}$ We can turn it into a diagonal matrix as follows, using the <code class=\"code-inline tex2jax_ignore\">diag<\/code> command from mpmath:<\/p><\/div>"]},
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":["S1 = mp.diag(S)\nS1"], "outputs":[]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-721\">In case you don't recognize the square root of 3 by its decimal approximation, we can check:<\/p><\/div>"]},
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":["S1*S1"], "outputs":[]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-722\">The matrix $\\Sigma_A$ is supposed to be $2\\times 3\\text{,}$ with an additional column of zeros. If we want to define this as a SymPy matrix, we can. Note that we are using the <code class=\"code-inline tex2jax_ignore\">diag<\/code> command from SymPy this time, rather than mpmath.<\/p><\/div>"]},
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":["S2 = sy.Matrix(S1)\nSigA = S2.row_join(sy.Matrix([0,0]))\nSigA"], "outputs":[]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-723\">We will see in a minute that this matrix is not really necessary, because the <code class=\"code-inline tex2jax_ignore\">svd<\/code> algorithm in mpmath is a bit different from the one in Nicholson.<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-724\">The matrix $P$ is already a SymPy matrix, as it turns out, but $Q$ is not.<\/p><\/div>"]},
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":["QT"], "outputs":[]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-725\">Notice that this matrix is not square! The output from the algorithm also has the transpose already applied.<\/p><\/div>"]},
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":["Q1 = QT.T\nQ1"], "outputs":[]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-726\">The matrix <code class=\"code-inline tex2jax_ignore\">Q1<\/code> we get here is (up to a difference in sign) the first two columns of the matrix $Q$ found in Nicholson. This makes sense: the matrix $S$ that we get from the algorithm is $D_A\\text{,}$ not $\\Sigma_A\\text{.}$ The third column of $\\Sigma_A$ is zero. So when we multiply by $Q^T\\text{,}$ the third row of $Q^T$ (that is, the third column of $Q$) is lost. That is, $S(QT)=\\Sigma_AQ^T\\text{.}$<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-727\">To confirm that everything worked, we can multiply:<\/p><\/div>"]},
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":["P*S1*QT"], "outputs":[]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-728\">Maybe that's not so enlightening, given all the decimal places. Let's check the difference with the matrix $A\\text{:}$<\/p><\/div>"]},
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":["A-P*S1*QT"], "outputs":[]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-729\">Looks like some pretty small numbers there. The mpmath library includes the <code class=\"code-inline tex2jax_ignore\">chop<\/code> command for truncating decimals.<\/p><\/div>"]},
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":["mp.chop(A-P*S1*QT)"], "outputs":[]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-730\">The Singular Value Decomposition has a lot of useful appplications, some of which are described in Nicholson's book. On a very fundamental level the <abbr xmlns:svg=\"http:\/\/www.w3.org\/2000\/svg\" class=\"initialism\">SVD<\/abbr> provides us with information on some of the most essential properties of the matrix $A\\text{,}$ and any system of equations with $A$ as its coefficient matrix.<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-731\">Recall the following definitions for an $m\\times n$ matrix $A\\text{:}$<\/p><ol class=\"decimal\"><li id=\"li-133\"><p id=\"p-732\">The <dfn xmlns:svg=\"http:\/\/www.w3.org\/2000\/svg\" class=\"terminology\">rank<\/dfn> of $A$ is the number of leadning ones in the <abbr class=\"initialism\">RREF<\/abbr> of $A\\text{,}$ which is also equal to the dimension of the column space of $A$ (or if you prefer, the dimension of $\\im (T_A)$).<\/p><\/li><li id=\"li-134\"><p id=\"p-733\">The <dfn class=\"terminology\">column space<\/dfn> of $A\\text{,}$ denoted $\\csp(A)\\text{,}$ is the subspace of $\\R^m$ spanned by the columns of $A\\text{.}$ (This is the image of the matrix transformation $T_A\\text{;}$ it is also the space of all vectors $\\mathbf{b}$ for which the system $A\\xx=\\mathbf{b}$ is consistent.)<\/p><\/li><li id=\"li-135\"><p id=\"p-734\">The <dfn class=\"terminology\">row space<\/dfn> of $A\\text{,}$ denoted $\\operatorname{row}(A)\\text{,}$ is the span of the rows of $A\\text{,}$ viewed as column vectors in $\\R^n\\text{.}$<\/p><\/li><li id=\"li-136\"><p id=\"p-735\">The <dfn class=\"terminology\">null space<\/dfn> of $A$ is the space of solutions to the homogeneous system $A\\xx=\\zer\\text{.}$ This is, of course, equal the kernel of the associated transformation $T_A\\text{.}$<\/p><\/li><\/ol><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-736\">There are some interesting relationships among these spaces, which are left as an exercise. (Solutions are in Nicholson if you get stumped.)<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><article class=\"exercise exercise-like\" id=\"exercise-60\"><h6 xmlns:svg=\"http:\/\/www.w3.org\/2000\/svg\" class=\"heading\"><span class=\"type\">Exercise<\/span><span class=\"space\"> <\/span><span class=\"codenumber\">4.5.3<\/span><span class=\"period\">.<\/span><\/h6><p id=\"p-737\">Let $A$ be an $m\\times n$ matrix. Prove the following:<\/p><ol class=\"decimal\"><li id=\"li-137\"><p id=\"p-738\">$\\displaystyle (\\operatorname{row}(A))^\\bot = \\nll(A)$<\/p><\/li><li id=\"li-138\"><p id=\"p-739\">$\\displaystyle (\\csp(A))^\\bot = \\nll(A^T)$<\/p><\/li><\/ol><\/article><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-740\">Here's the cool thing about the <abbr xmlns:svg=\"http:\/\/www.w3.org\/2000\/svg\" class=\"initialism\">SVD<\/abbr>. Let $\\sigma_1\\geq \\sigma_2\\geq \\cdots \\geq \\sigma_r\\gt 0$ be the positive singular values of $A\\text{.}$ Let $\\vecq_1,\\ldots, \\vecq_r,\\ldots, \\vecq_n$ be the orthonormal basis of eigenvectors for $A^TA\\text{,}$ and let $\\vecp_1,\\ldots, \\vecp_r,\\ldots, \\vecp_m$ be the orthonormal basis of $\\R^m$ constructed in the <abbr class=\"initialism\">SVD<\/abbr> algorithm. Then:<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><ol id=\"p-741\" class=\"decimal\"><li id=\"li-139\"><p id=\"p-742\">$\\displaystyle \\rank(A)=r$<\/p><\/li><li id=\"li-140\"><p id=\"p-743\">$\\vecq_1,\\ldots, \\vecq_r$ form a basis for $\\operatorname{row}(A)\\text{.}$<\/p><\/li><li id=\"li-141\"><p id=\"p-744\">$\\vecp_1,\\ldots, \\vecp_r$ form a basis for $\\csp(A)$ (and thus, the “row rank” and “column rank” of $A$ are the same).<\/p><\/li><li id=\"li-142\"><p id=\"p-745\">$\\vecq_{r+1},\\ldots, \\vecq_n$ form a basis for $\\nll(A)\\text{.}$ (And these are therefore the basis solutions of $A\\xx=\\zer\\text{!}$)<\/p><\/li><li id=\"li-143\"><p id=\"p-746\">$\\vecp_{r+1},\\ldots, \\vecp_m$ form a basis for $\\nll(A^T)\\text{.}$<\/p><\/li><\/ol><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-747\">If you want to explore this further, have a look at the excellent <a class=\"external\" href=\"https:\/\/www.juanklopper.com\/wp-content\/uploads\/2015\/03\/III%5F05%5FSingular%5Fvalue%5Fdecomposition.html\" target=\"_blank\">notebook by Dr. Juan H Klopper<\/a>. The <code class=\"code-inline tex2jax_ignore\">ipynb<\/code> file can be found <a class=\"external\" href=\"https:\/\/github.com\/juanklopper\/MIT%5FOCW%5FLinear%5FAlgebra%5F18%5F06\" target=\"_blank\">on his GitHub page<\/a>. In it, he takes you through various approaches to finding the singular value decomposition, using the method above, as well as using NumPy and SciPy (which, for industrial applications, are superior to SymPy and mpmath).<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><h6 class=\"heading\"><span class=\"title\">Polar Decomposition.<\/span><\/h6><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-748\">For any $n\\times n$ matrix $A\\text{,}$ there exists an orthogonal (or unitary) matrix $P$ such that<\/p><div xmlns:svg=\"http:\/\/www.w3.org\/2000\/svg\" class=\"displaymath\">\n\\begin{equation*}\nA = P\\sqrt{A^TA}\\text{.}\n\\end{equation*}\n<\/div><p class=\"continuation\">This is meant to remind you of the polar decomposition<\/p><div class=\"displaymath\">\n\\begin{equation*}\nz = e^{i\\theta}\\sqrt{\\bar{z}z}\n\\end{equation*}\n<\/div><p class=\"continuation\">for a complex number.<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-749\">One way to compute the polar decomposition is using the Singular Value Decomposition (see Nicholson's text). Note that both $P$ and $\\sqrt{A^TA}$ can be diagonalized, but usually not with the same orthonormal basis.<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><h6 class=\"heading\"><span class=\"title\">QR Factorization.<\/span><\/h6><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-750\">Suppose $A$ is an $m\\times n$ matrix with independent columns. (Question: for this to happen, which is true — $m\\geq n\\text{,}$ or $n\\geq m\\text{?}$)<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-751\">A $QR$-factorization of $A$ is a factorization of the form $A=QR\\text{,}$ where $Q$ is $m\\times n\\text{,}$ with orthonormal columns, and $R$ is an invertible upper-triangular ($n\\times n$) matrix with positive diagonal entries. If $A$ is a square matrix, $Q$ will be orthogonal.<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-752\">A lot of the methods we're looking at here involve more sophisticated numerical techniques than SymPy is designed to handle. If we wanted to spend time on these topics, we'd have to learn a bit about the NumPy package, which has built in tools for finding things like polar decomposition and singular value decomposition. However, SymPy does know how to do $QR$ factorization. After defining a matrix <code class=\"code-inline tex2jax_ignore\">A<\/code>, we can use the command<\/p><pre class=\"code-display tex2jax_ignore\">\n          Q, R = A.QRdecomposition()\n        <\/pre><p class=\"continuation\">.<\/p><\/div>"]},
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":["from sympy import Matrix,init_printing\ninit_printing()\nA = Matrix(3,3,[1,-2,3,3,-1,2,4,2,5])\nQ, R = A.QRdecomposition()\nA, Q, R"], "outputs":[]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-753\">Let's check that the matrix $Q$ really is orthogonal:<\/p><\/div>"]},
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":["Q**(-1) == Q.T"], "outputs":[]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-754\">Details of how to perform the QR factorization can be found in Nicholson's textbook. It's essentially a consequence of performing the Gram-Schmidt algorithm on the columns of $A\\text{,}$ and keeping track of our work.<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-755\">The calculation above is a symbolic computation, which is nice for understanding what's going on. The reason why the $QR$ factorization is useful in practice is that there are efficient numerical methods for doing it (with good control over rounding errors). Our next topic looks at a useful application of the $QR$ factorization.<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><h6 class=\"heading hide-type\"><span xmlns:svg=\"http:\/\/www.w3.org\/2000\/svg\" class=\"type\">Subsection<\/span> <span class=\"codenumber\">4.5.2<\/span> <span class=\"title\">Computing Eigenvalues<\/span><\/h6><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-756\">Our first method focuses on the dominant eigenvalue of a matrix. An eigenvalue is dominant if it is larger in absolute value than all other eigenvalues. For example, if $A$ has eigenvalues $1,3,-2,-5\\text{,}$ then $-5$ is the dominant eigenvalue.<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-757\">If $A$ has eigenvalues $1,3,0,-4,4$ then there is no dominant eigenvalue. Any eigenvector corresponding to a dominant eigenvalue is called a dominant eigenvector.<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><h6 class=\"heading\"><span class=\"title\">The Power Method.<\/span><\/h6><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-758\">If a matrix $A$ has a dominant eigenvalue, there is a method for finding it (approximately) that does not involve finding and factoring the characteristic polynomial of $A\\text{.}$<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-759\">We start with some initial guess $x_0$ for a dominant eigenvector. We then set $x_{k+1} = Ax_k$ for each $k\\geq 0\\text{,}$ giving a sequence<\/p><div xmlns:svg=\"http:\/\/www.w3.org\/2000\/svg\" class=\"displaymath\">\n\\begin{equation*}\nx_0, Ax_0, A^2x_0, A^3x_0,\\ldots\\text{.}\n\\end{equation*}\n<\/div><p class=\"continuation\">We expect (for reasons we'll explain) that $\\lVert x_k-x\\rVert \\to 0$ as $k\\to\\infty\\text{,}$ where $x$ is a dominant eigenvector. Let's try an example.<\/p><\/div>"]},
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":["A = Matrix(2,2,[1,-4,-3,5])\nA,A.eigenvects()"], "outputs":[]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-760\">The dominant eigenvalue is $\\lambda = 7\\text{.}$ Let's try an initial guess of $x_0=\\begin{bmatrix}1\\\\0\\end{bmatrix}$ and see what happens.<\/p><\/div>"]},
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":["x0 = Matrix(2,1,[1,0])\nL = list()\nfor k in range(10):\n    L.append(A**k*x0)\nL"], "outputs":[]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-761\">We might want to confirm whether that rather large fraction is close to $\\frac23\\text{.}$ To do so, we can get the computer to divide the numerator by the denominator.<\/p><\/div>"]},
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":["L[9][0]\/L[9][1]"], "outputs":[]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-762\">The above might show you the fraction rather than its decimal approximation. (This may depend on whether you're on Sage or Jupyter.) To get the decimal, try wrapping the above in <code class=\"code-inline tex2jax_ignore\">float()<\/code> (or <code class=\"code-inline tex2jax_ignore\">N<\/code>, or append with <code class=\"code-inline tex2jax_ignore\">.evalf()<\/code>).<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-763\">For the eigenvalue, we note that if $Ax=\\lambda x\\text{,}$ then<\/p><div xmlns:svg=\"http:\/\/www.w3.org\/2000\/svg\" class=\"displaymath\">\n\\begin{equation*}\n\\frac{x\\cdot Ax}{\\lVert x\\rVert^2} = \\frac{x\\cdot (\\lambda x)}{\\lVert x\\rVert^2} = \\lambda\\text{.}\n\\end{equation*}\n<\/div><p class=\"continuation\">This leads us to consider the Rayleigh quotients<\/p><div class=\"displaymath\">\n\\begin{equation*}\nr_k = \\frac{x_k\\cdot x_{k+1}}{\\lVert x_k\\rVert^2}\\text{.}\n\\end{equation*}\n<\/div><\/div>"]},
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":["M = list()\nfor k in range(9):\n    M.append((L[k].dot(L[k+1]))\/(L[k].dot(L[k])))\nM"], "outputs":[]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-764\">We can convert a rational number r to a float using either <code class=\"code-inline tex2jax_ignore\">N(r)<\/code> or <code class=\"code-inline tex2jax_ignore\">r.evalf()<\/code>. (The latter seems to be the better bet when working with a list.)<\/p><\/div>"]},
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":["M2 = list()\nfor k in range(9):\n    M2.append((M[k]).evalf())\nM2"], "outputs":[]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><h6 class=\"heading\"><span class=\"title\">The QR Algorithm.<\/span><\/h6><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-765\">Given an $n\\times n$ matrix $A\\text{,}$ we know we can write $A=QR\\text{,}$ with $Q$ orthogonal and $R$ upper-triangular. The $QR$-algorithm exploits this fact. We set $A_1=A\\text{,}$ and write $A_1=Q_1R_1\\text{.}$<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-766\">Then we set $A_2 = R_1Q_1\\text{,}$ and factor: $A_2=Q_2R_2\\text{.}$ Notice $A_2 = R_1Q_1 = Q_1^TA_1Q_1\\text{.}$ Since $A_2$ is similar to $A_1\\text{,}$ $A_2$ has the same eigenvalues as $A_1=A\\text{.}$<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-767\">Next, set $A_3 = R_2Q_2\\text{,}$ and factor as $A_3 = Q_3R_3\\text{.}$ Since $A_3 = Q_2^TA_2Q_2\\text{,}$ $A_3$ has the same eigenvalues as $A_2\\text{.}$ In fact, $A_3 = Q_2^T(Q_1^TAQ_1)Q_2 = (Q_1Q_2)^TA(Q_1Q_2)\\text{.}$<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-768\">After $k$ steps we have $A_{k+1} = (Q_1\\cdots Q_k)^TA(Q_1\\cdots Q_k)\\text{,}$ which still has the same eigenvalues as $A\\text{.}$ By some sort of dark magic, this sequence of matrices converges to an upper triangular matrix with eigenvalues on the diagonal!<\/p><\/div>"]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-769\">Consider the matrix $A = \\begin{bmatrix}5&-2&3\\\\0&4&0\\\\0&-1&3\\end{bmatrix}$<\/p><\/div>"]},
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":["A = Matrix(3,3,[5,-2,3,0,4,0,0,-1,3])\nA.eigenvals()"], "outputs":[]},
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":["Q1,R1 = A.QRdecomposition()\nA2=R1*Q1\nA2,Q1,R1"], "outputs":[]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-770\">Now we repeat the process:<\/p><\/div>"]},
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":["Q2,R2 = A2.QRdecomposition()\nA3=R2*Q2\nA3.evalf()"], "outputs":[]},
{"cell_type":"markdown", "metadata":{}, "source":["<div class=\"mathbook-content\"><p id=\"p-771\">Do this a few more times, and see what results! (If someone can come up with a way to code this as a loop, let me know!) The diagonal entries should get closer to $5,4,3\\text{,}$ respectively, and the $(3,2)$ entry should get closer to $0\\text{.}$<\/p><\/div>"]},
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":[""], "outputs":[]},
{"cell_type":    "code", "execution_count":null, "metadata":{}, "source":[""], "outputs":[]}
],
"nbformat": 4, "nbformat_minor": 0, "metadata": {"kernelspec": {"display_name": "", "name": "python3"}, "language_info": {"codemirror_mode": {"name": "ipython", "version": 3}, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.4"}, "name": "section-matrix-factor.ipynb"}
}