📚 The CoCalc Library - books, templates and other resources
License: OTHER
%auto %html(hide=True) <div class="mathbook-content"> <link rel="stylesheet" type="text/css" href="http://buzzard.ups.edu/mathbook-content.css"> <link rel="stylesheet" type="text/css" href="https://aimath.org/mathbook/mathbook-add-on.css"> </div>
%auto %html(hide=True) <div class="mathbook-content"><div class="hidden-content" style="display:none">\( \newcommand{\lt}{<} \newcommand{\gt}{>} \newcommand{\amp}{&} \)</div></div>
%auto %html(hide=True) <div class="mathbook-content"><table width="90%" style="font-size: 200%;"><tr></tr></table></div>
%auto %html(hide=True) <div class="mathbook-content"><section class="article" id="PDM"></section></div>
%auto %html(hide=True) <div class="mathbook-content"><section class="frontmatter" id="frontmatter-1"></section></div>
%auto %html(hide=True) <div class="mathbook-content"> <h2 class="heading"> <span class="title">Sage and Linear Algebra Worksheet:</span> <span class="subtitle">FCLA Section PDM</span> </h2> <div class="author"> <div class="author-name">Robert Beezer</div> <div class="author-info">Department of Mathematics and Computer Science<br>University of Puget Sound</div> </div> <div class="date">Fall 2019</div> </div>
%auto %html(hide=True) <div class="mathbook-content"><section class="section" id="section-1"><h6 class="heading hide-type"> <span class="type">Section</span> <span class="codenumber">1</span> <span class="title">LU Decomposition, Triangular Form</span> </h6></section></div>
%auto %html(hide=True) <div class="mathbook-content"><p id="p-1">This is a topic not covered in our text. You <em class="emphasis">can</em> find a discussion in <span class="booktitle">A Second Course in Linear Algebra</span> at <a class="external" href="http://linear.ups.edu/scla/html/index.html" target="_blank"><code class="code-inline tex2jax_ignore">http://linear.ups.edu/scla/html/index.html</code></a>.</p></div>
%auto %html(hide=True) <div class="mathbook-content"> <p id="p-2">Our goal is to row-reduce a matrix with elementary matrices, track the changes, and arrive at an expression for a square matrix \(A\) as a product of a lower-triangular matrix, \(L\text{,}\) and an upper-triangular matrix, \(U\text{,}\) that is</p> <div class="displaymath"> \begin{equation*} A=LU \end{equation*} </div> <p>the so-called <dfn class="terminology">LU decomposition</dfn>. I sometimes prefer to call it <dfn class="terminology">triangular form</dfn>.</p> </div>
%auto %html(hide=True) <div class="mathbook-content"><p id="p-3">There are no exercises in this worksheet, but instead there is a careful and detailed exposition of using elementary matrices (row operations) to arrive at a <dfn class="terminology">matrix decomposition</dfn>. There are many kinds of matrix decompositions, such as the <dfn class="terminology">singular value decomposition</dfn> (SVD). Five or six such decompositions form a central part of the linear algebra canon. Again, see <span class="booktitle">A Second Course in Linear Algebra</span> for details on these.</p></div>
%auto %html(hide=True) <div class="mathbook-content"><p id="p-4">We decompose a \(5\times 5\) matrix. It is most natural to describe an LU decomposition of a square matrix, but the decomposition can be generalized to rectangular matrices.</p></div>
A = matrix(QQ, [[-6, -10, 0, 10, 14], [ 2, 3, 0, -4, -3], [ 0, -2, -3, 1, 8], [ 5, 6, -3, -7, -3], [-1, 1, 6, -1, -8]]) A
%auto %html(hide=True) <div class="mathbook-content"><p id="p-5">Elementary matrices to “do” row operations in the first column.</p></div>
actionA = elementary_matrix(QQ, 5, row1=1, row2=0, scale=-2)*elementary_matrix(QQ, 5, row1=3, row2=0, scale=-5)*elementary_matrix(QQ, 5, row1=4, row2=0, scale=1)*elementary_matrix(QQ, 5, row1=0, scale=-1/6) B = actionA*A B
%auto %html(hide=True) <div class="mathbook-content"><p id="p-6">Now in second column, moving to <dfn class="terminology">row-echelon form</dfn> (i.e. <em class="emphasis">not</em> <dfn class="terminology">reduced row-echelon form</dfn>).</p></div>
actionB = elementary_matrix(QQ, 5, row1=2, row2=1, scale=2)*elementary_matrix(QQ, 5, row1=3, row2=1, scale=7/3)*elementary_matrix(QQ, 5, row1=4, row2=1, scale=-8/3)*elementary_matrix(QQ, 5, row1=1, scale=-3) C = actionB*B C
%auto %html(hide=True) <div class="mathbook-content"><p id="p-7">The “bottom” of the third column.</p></div>
actionC = elementary_matrix(QQ, 5, row1=3, row2=2, scale=3)*elementary_matrix(QQ, 5, row1=4, row2=2, scale=-6)*elementary_matrix(QQ, 5, row1=2, scale=-1/3) D = actionC*C D
%auto %html(hide=True) <div class="mathbook-content"><p id="p-8">And now the penultimate column.</p></div>
actionD = elementary_matrix(QQ, 5, row1=4, row2=3, scale=-2)*elementary_matrix(QQ, 5, row1=3, scale=1) E = actionD*D E
%auto %html(hide=True) <div class="mathbook-content"><p id="p-9">And done.</p></div>
actionE = elementary_matrix(QQ, 5, row1=4, scale=1) F = actionE*E F
%auto %html(hide=True) <div class="mathbook-content"> <p id="p-10">Clearly, <code class="code-inline tex2jax_ignore">F</code> has determinant 1, since it is an upper triangular matrix with diagonal entries equal to \(1\text{.}\) By tracking the effect of the above manipulations (tantamount to performing row operations) we expect that</p> <div class="displaymath"> \begin{equation*} \det(A) = \left(\frac{1}{-1/6}\right)\left(\frac{1}{-3}\right)\left(\frac{1}{-1/3}\right)\left(\frac{1}{1}\right)\left(\frac{1}{1}\right)\det(F) = -6. \end{equation*} </div> <p>Let's check.</p> </div>
A.determinant()
%auto %html(hide=True) <div class="mathbook-content"><p id="p-11">Yep. But it gets better. <code class="code-inline tex2jax_ignore">F</code> is the product of the “action” matrices on the left of <code class="code-inline tex2jax_ignore">A</code>.</p></div>
total_action = prod([actionE, actionD, actionC, actionB, actionA]) total_action
%auto %html(hide=True) <div class="mathbook-content"><p id="p-12">Notice that the elementary matrices we used are all lower triangular (because we just formed zeros below the diagonal of the original matrix as we brought it to row-echelon form, and there were no row swaps). Hence their product is again lower triangular. Now check that we have the correct matrix.</p></div>
F == total_action * A
%auto %html(hide=True) <div class="mathbook-content"><p id="p-13">The “total action” matrix is a product of elementary matrices, which are individually nonsingular. So their product is nonsingular. Futhermore, the inverse is again lower triangular.</p></div>
ta_inv = total_action.inverse() ta_inv
%auto %html(hide=True) <div class="mathbook-content"><p id="p-14">We reach our goal by rearranging the equality above, writing <code class="code-inline tex2jax_ignore">A</code> as a product of a lower-triangular matrix with an upper-triangular matrix.</p></div>
A == ta_inv * F
%auto %html(hide=True) <div class="mathbook-content"><p id="p-15">Yes! So we have decomposed the original matrix (<code class="code-inline tex2jax_ignore">A</code>) into the product of a lower triangular matrix (inverse of the total action matrix) and an upper triangular matrix with all ones on the diagonal (<code class="code-inline tex2jax_ignore">F</code>, the original matrix in row-echelon form).</p></div>
A, ta_inv, F
%auto %html(hide=True) <div class="mathbook-content"><p id="p-16">This decomposition (the <dfn class="terminology">LU decomposition</dfn>) can be useful for solving systems quickly. You <dfn class="terminology">forward solve</dfn> with \(L\text{,}\) then <dfn class="terminology">back solve</dfn> with \(U\text{.}\)</p></div>
%auto %html(hide=True) <div class="mathbook-content"><p id="p-17">More specifically, suppose you want to solve \(A\mathbf{x}=\mathbf{b}\) for \(\mathbf{x}\text{,}\) and you have a decomposition \(A=LU\text{.}\) First solve the intermediate system, \(L\mathbf{y}=\mathbf{b}\) for \(\mathbf{y}\text{,}\) which can be accomplished easily by determining the entries of \(\mathbf{y}\) in order, exploiting the lower triangular nature of \(L\text{.}\) This is what is meant by the term <dfn class="terminology">forward solve</dfn>.</p></div>
%auto %html(hide=True) <div class="mathbook-content"><p id="p-18">With a solution for \(\mathbf{y}\text{,}\) form the system \(U\mathbf{x}=\mathbf{y}\text{.}\) You can check that a solution, \(\mathbf{x}\text{,}\) to this system is also a solution to the original system \(A\mathbf{x}=\mathbf{b}\text{.}\) Further, this solution can be found easily by determining the entries of \(\mathbf{x}\) in reverse order, exploiting the upper triangular nature of \(U\text{.}\) This is what is meant by the term <dfn class="terminology">back solve</dfn>.</p></div>
%auto %html(hide=True) <div class="mathbook-content"><p id="p-19">We solve <em class="emphasis">two</em> simple systems, but only do half as many row-operations as if we went fully to reduced row-echelon form. If you count the opertions carefully, you will see that this is a big win, roughly reducing computation time by a factor of half for large systems.</p></div>
%auto %html(hide=True) <div class="mathbook-content"><article class="conclusion" id="conclusion-1"><h5 class="heading"><span></span></h5></article></div>
%auto %html(hide=True) <div class="mathbook-content"><p id="p-20">This work is Copyright 2016–2019 by Robert A. Beezer. It is licensed under a <a class="external" href="https://creativecommons.org/licenses/by-sa/4.0/" target="_blank">Creative Commons Attribution-ShareAlike 4.0 International License</a>.</p></div>
%auto %html(hide=True) <div class="mathbook-content"><table width="90%" style="font-size: 200%;"><tr></tr></table></div>