**Authors**
Robert A. Beezer

**License**
GFDL-1.2-no-invariants-or-later

Sage Primer for Linear Algebra A First Course in Linear Algebra Robert A. Beezer University of Puget Sound Version 3.50 (December 30, 2015) c 20042015 Robert A. Beezer Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation; with no Invariant Sections, no Front- Cover Texts, and no Back-Cover Texts. A copy of the license is included in the appendix entitled “GNU Free Documentation License”. Contents iii Chapter SLE Systems of Linear Equations Section SSLE: Solving Systems of Linear Equations GS Getting Started Sage is a powerful system for studying and exploring many different areas of math- ematics. In the next section, and the majority of the remaining section, we will include short descriptions and examples using Sage. You can read a bit more about Sage in the Preface. If you are not already reading this in an electronic version, you may want to investigate obtaining the worksheet version of this book, where the examples are “live” and editable. Most of your interaction with Sage will be by typing commands into a compute cell. That is a compute cell just below this paragraph. Click once inside the compute cell and you will get a more distinctive border around it, a blinking cursor inside, plus a cute little “evaluate” link below it. At the cursor, type 2+2 and then click on the evaluate link. Did a 4 appear below the cell? If so, you have successfully sent a command off for Sage to evaluate and you have received back the (correct) answer. Here is another compute cell. Try evaluating the command factorial(300). Hmmmmm. That is quite a big integer! The slashes you see at the end of each line mean the result is continued onto the next line, since there are 615 digits in the result. To make new compute cells, hover your mouse just above another compute cell, or just below some output from a compute cell. When you see a skinny blue bar across the width of your worksheet, click and you will open up a new compute cell, ready for input. Note that your worksheet will remember any calculations you make, in the order you make them, no matter where you put the cells, so it is best to stay organized and add new cells at the bottom. Try placing your cursor just below the monstrous value of 300! that you have. Click on the blue bar and try another factorial computation in the new compute cell. Each compute cell will show output due to only the very last command in the cell. Try to predict the following output before evaluating the cell. a = 10 b = 6 a = a + 20 a 30 The following compute cell will not print anything since the one command does not create output. But it will have an effect, as you can see when you execute 1 Chapter SLE 2 the subsequent cell. Notice how this uses the value of b from above. Execute this compute cell once. Exactly once. Even if it appears to do nothing. If you execute the cell twice, your credit card may be charged twice. b = b + 50 Now execute this cell, which will produce some output. b + 20 76 So b came into existence as 6. Then a cell added 50. This assumes you only executed this cell once! In the last cell we create b+20 (but do not save it) and it is this value that is output. You can combine several commands on one line with a semi-colon. This is a great way to get multiple outputs from a compute cell. The syntax for building a matrix should be somewhat obvious when you see the output, but if not, it is not particularly important to understand now. f(x) = x^8 - 7*x^4; f x |--> x^8 - 7*x^4 f; print ; f.derivative() x |--> x^8 - 7*x^4 <BLANKLINE> x |--> 8*x^7 - 28*x^3 g = f.derivative() g.factor() 4*(2*x^4 - 7)*x^3 Some commands in Sage are “functions,” an example is factorial() above. Other commands are “methods” of an object and are like characteristics of objects, examples are .factor() and .derivative() as methods of a function. To comment on your work, you can open up a small word-processor. Hover your mouse until you get the skinny blue bar again, but now when you click, also hold the SHIFT key at the same time. Experiment with fonts, colors, bullet lists, etc and then click the “Save changes” button to exit. Double-click on your text if you need to go back and edit it later. Open the word-processor again to create a new bit of text (maybe next to the empty compute cell just below). Type all of the following exactly, but do not include any backslashes that might precede the dollar signs in the print version: Pythagorean Theorem: \$c^2=a^2+b^2\$ and save your changes. The symbols between the dollar signs are written accord- ing to the mathematical typesetting language known as TeX — cruise the internet to learn more about this very popular tool. (Well, it is extremely popular among mathematicians and physical scientists.) Much of our interaction with sets will be through Sage lists. These are not really sets — they allow duplicates, and order matters. But they are so close to sets, and so easy and powerful to use that we will use them regularly. We will use a fun made-up list for practice, the quote marks mean the items are just text, with no special mathematical meaning. Execute these compute cells as we work through them. Chapter SLE 3 zoo = [’snake’, ’parrot’, ’elephant’, ’baboon’, ’beetle’] zoo [’snake’, ’parrot’, ’elephant’, ’baboon’, ’beetle’] So the square brackets define the boundaries of our list, commas separate items, and we can give the list a name. To work with just one element of the list, we use the name and a pair of brackets with an index. Notice that lists have indices that begin counting at zero. This will seem odd at first and will seem very natural later. zoo[2] ’elephant’ We can add a new creature to the zoo, it is joined up at the far right end. zoo.append(’ostrich’); zoo [’snake’, ’parrot’, ’elephant’, ’baboon’, ’beetle’, ’ostrich’] We can remove a creature. zoo.remove(’parrot’) zoo [’snake’, ’elephant’, ’baboon’, ’beetle’, ’ostrich’] We can extract a sublist. Here we start with element 1 (the elephant) and go all the way up to, but not including, element 3 (the beetle). Again a bit odd, but it will feel natural later. For now, notice that we are extracting two elements of the lists, exactly 3 − 1 = 2 elements. mammals = zoo[1:3] mammals [’elephant’, ’baboon’] Often we will want to see if two lists are equal. To do that we will need to sort a list first. A function creates a new, sorted list, leaving the original alone. So we need to save the new one with a new name. newzoo = sorted(zoo) newzoo [’baboon’, ’beetle’, ’elephant’, ’ostrich’, ’snake’] zoo.sort() zoo [’baboon’, ’beetle’, ’elephant’, ’ostrich’, ’snake’] Notice that if you run this last compute cell your zoo has changed and some commands above will not necessarily execute the same way. If you want to experi- ment, go all the way back to the first creation of the zoo and start executing cells again from there with a fresh zoo. A construction called a “list comprehension” is especially powerful, especially since it almost exactly mirrors notation we use to describe sets. Suppose we want to form the plural of the names of the creatures in our zoo. We build a new list, based on all of the elements of our old list. Chapter SLE 4 plurality_zoo = [animal+’s’ for animal in zoo] plurality_zoo [’baboons’, ’beetles’, ’elephants’, ’ostrichs’, ’snakes’] Almost like it says: we add an “s” to each animal name, for each animal in the zoo, and place them in a new list. Perfect. (Except for getting the plural of “ostrich” wrong.) One final type of list, with numbers this time. The range() function will create lists of integers. In its simplest form an invocation like range(12) will create a list of 12 integers, starting at zero and working up to, but not including, 12. Does this sound familiar? dozen = range(12); dozen [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11] Here are two other forms, that you should be able to understand by studying the examples. teens = range(13, 20); teens [13, 14, 15, 16, 17, 18, 19] decades = range(1900, 2000, 10); decades [1900, 1910, 1920, 1930, 1940, 1950, 1960, 1970, 1980, 1990] There is a “Save” button in the upper-right corner of your worksheet. This will save a current copy of your worksheet that you can retrieve from within your notebook again later, though you have to re-execute all the cells when you re-open the worksheet later. There is also a “File” drop-down list, on the left, just above your very top compute cell (not be confused with your browser’s File menu item!). You will see a choice here labeled “Save worksheet to a file...” When you do this, you are creating a copy of your worksheet in the “sws” format (short for “Sage WorkSheet”). You can email this file, or post it on a website, for other Sage users and they can use the “Upload” link on their main notebook page to incorporate a copy of your worksheet into their notebook. There are other ways to share worksheets that you can experiment with, but this gives you one way to share any worksheet with anybody almost anywhere. We have covered a lot here in this section, so come back later to pick up tidbits you might have missed. There are also many more features in the notebook that we have not covered. Section RREF: Reduced Row-Echelon Form M Matrices Matrices are fundamental objects in linear algebra and in Sage, so there are a variety of ways to construct a matrix in Sage. Generally, you need to specify what types of entries the matrix contains (more on that in a minute), the number of rows and columns, and the entries themselves. First, let us dissect an example: Chapter SLE 5 A = matrix(QQ, 2, 3, [[1, 2, 3], [4, 5, 6]]) A [1 2 3] [4 5 6] QQ is the set of all rational numbers (fractions with an integer numerator and denominator), 2 is the number of rows, 3 is the number of columns. Sage under- stands a list of items as delimited by brackets ([,]) and the items in the list can again be lists themselves. So [[1, 2, 3], [4, 5, 6]] is a list of lists, and in this context the inner lists are rows of the matrix. There are various shortcuts you can employ when creating a matrix. For exam- ple, Sage is able to infer the size of the matrix from the lists of entries. B = matrix(QQ, [[1, 2, 3], [4, 5, 6]]) B [1 2 3] [4 5 6] Or you can specify how many rows the matrix will have and provide one big grand list of entries, which will get chopped up, row by row, if you prefer. C = matrix(QQ, 2, [1, 2, 3, 4, 5, 6]) C [1 2 3] [4 5 6] It is possible to also skip specifying the type of numbers used for entries of a matrix, however this is fraught with peril, as Sage will make an informed guess about your intent. Is this what you want? Consider when you enter the single character “2” into a computer program like Sage. Is this the integer 2, the rational number 21 , the real number 2.00000, the complex number 2 + 0i, or the polynomial p(x) = 2? In context, us humans can usually figure it out, but a literal-minded computer is not so smart. It happens that the operations we can perform, and how they behave, are influenced by the type of the entries in a matrix. So it is important to get this right and our advice is to be explicit and be in the habit of always specifying the type of the entries of a matrix you create. Mathematical objects in Sage often come from sets of similar objects. This set is called the “parent” of the element. We can use this to learn how Sage deduces the type of entries in a matrix. Execute the following three compute cells in the Sage notebook, and notice how the three matrices are constructed to have entries from the integers, the rationals and the reals. A = matrix(2, 3, [[1, 2, 3], [4, 5, 6]]) A.parent() Full MatrixSpace of 2 by 3 dense matrices over Integer Ring B = matrix(2, 3, [[1, 2/3, 3], [4, 5, 6]]) B.parent() Chapter SLE 6 Full MatrixSpace of 2 by 3 dense matrices over Rational Field C = matrix(2, 3, [[1, sin(2.2), 3], [4, 5, 6]]) C.parent() Full MatrixSpace of 2 by 3 dense matrices over Real Field with 53 bits of precision Sage knows a wide variety of sets of numbers. These are known as “rings” or “fields” (see Section F), but we will call them “number systems” here. Examples include: ZZ is the integers, QQ is the rationals, RR is the real numbers with reason- able precision, and CC is the complex numbers with reasonable precision. We will present the theory of linear algebra over the complex numbers. However, in any computer system, there will always be complications surrounding the inability of digital arithmetic to accurately represent all complex numbers. In contrast, Sage can represent rational numbers exactly as the quotient of two (perhaps very large) integers. So our Sage examples will begin by using QQ as our number system and we can concentrate on understanding the key concepts. Once we have constructed a matrix, we can learn a lot about it (such as its parent). Sage is largely object-oriented, which means many commands apply to an object by using the “dot” notation. A.parent() is an example of this syntax, while the constructor matrix([[1, 2, 3], [4, 5, 6]]) is an exception. Here are a few examples, followed by some explanation: A = matrix(QQ, 2, 3, [[1,2,3],[4,5,6]]) A.nrows(), A.ncols() (2, 3) A.base_ring() Rational Field A[1,1] 5 A[1,2] 6 The number of rows and the number of columns should be apparent, .base_ring() gives the number system for the entries, as included in the information provided by .parent(). Computer scientists and computer languages prefer to begin counting from zero, while mathematicians and written mathematics prefer to begin counting at one. Sage and this text are no exception. It takes some getting used to, but the reasons for counting from zero in computer programs soon becomes very obvious. Counting from one in mathematics is historical, and unlikely to change anytime soon. So above, the two rows of A are numbered 0 and 1, while the columns are numbered 0, 1 and 2. So A[1,2] refers to the entry of A in the second row and the third column, i.e. 6. There is much more to say about how Sage works with matrices, but this is already a lot to digest. Use the space below to create some matrices (different ways) and examine them and their properties (size, entries, number system, parent). Chapter SLE 7 V Vectors Vectors in Sage are built, manipulated and interrogated in much the same way as matrices (see Sage M). However as simple lists (“one-dimensional,” not “two- dimensional” such as matrices that look more tabular) they are simpler to construct and manipulate. Sage will print a vector across the screen, even if we wish to interpret it as a column vector. It will be delimited by parentheses ((,)) which allows us to distinguish a vector from a matrix with just one row, if we look carefully. The number of “slots” in a vector is not referred to in Sage as rows or columns, but rather by “degree.” Here are some examples (remember to start counting from zero): v = vector(QQ, 4, [1, 1/2, 1/3, 1/4]) v (1, 1/2, 1/3, 1/4) v.degree() 4 v.parent() Vector space of dimension 4 over Rational Field v[2] 1/3 w = vector([1, 2, 3, 4, 5, 6]) w (1, 2, 3, 4, 5, 6) w.degree() 6 w.parent() Ambient free module of rank 6 over the principal ideal domain Integer Ring w[3] Chapter SLE 8 4 Notice that if you use commands like .parent() you will sometimes see refer- ences to “free modules.” This is a technical generalization of the notion of a vector, which is beyond the scope of this course, so just mentally convert to vectors when you see this term. The zero vector is super easy to build, but be sure to specify what number system your zero is from. z = zero_vector(QQ, 5) z (0, 0, 0, 0, 0) Notice that while using Sage, we try to remain consistent with our mathematical notation conventions. This is not required, as you can give items in Sage very long names if you wish. For example, the following is perfectly legitimate, as you can see. blatzo = matrix(QQ, 2, [1, 2, 3, 4]) blatzo [1 2] [3 4] In fact, our use of capital letters for matrices actually contradicts some of the conventions for naming objects in Sage, so there are good reasons for not mirroring our mathematical notation. AM Augmented Matrix Sage has a matrix method, .augment(), that will join two matrices, side-by-side provided they both have the same number of rows. The same method will allow you to augment a matrix with a column vector, as described in Definition AM, provided the number of entries in the vector matches the number of rows for the matrix. Here we reprise the construction in Example AMAA. We will now format our matrices as input across several lines, a practice you may use in your own worksheets, or not. A = matrix(QQ, 3, 3, [[1, -1, 2], [2, 1, 1], [1, 1, 0]]) b = vector(QQ, [1, 8, 5]) M = A.augment(b) M [ 1 -1 2 1] [ 2 1 1 8] [ 1 1 0 5] Notice that the matrix method .augment() needs some input, in the above case, the vector b. This will explain the need for the parentheses on the end of the “dot” commands, even if the particular command does not expect input. Some methods allow optional input, typically using keywords. Matrices can track subdivisions, making breaks between rows and/or columns. When augment- ing, you can ask for the subdivision to be included. Evalute the compute cell above if you have not already, so that A and b are defined, and then evaluate: Chapter SLE 9 M = A.augment(b, subdivide=True) M [ 1 -1 2| 1] [ 2 1 1| 8] [ 1 1 0| 5] As a partial demonstration of manipulating subdivisions of matrices we can reset the subdivisions of M with the .subdivide() method. We provide a list of rows to subdivide before, then a list of columns to subdivide before, where we remember that counting begins at zero. M.subdivide([1,2],[1]) M [ 1|-1 2 1] [--+--------] [ 2| 1 1 8] [--+--------] [ 1| 1 0 5] RO Row Operations Sage will perform individual row operations on a matrix. This can get a bit tedious, but it is better than doing the computations (wrong, perhaps) by hand, and it can be useful when building up more complicated procedures for a matrix. For each row operation, there are two similar methods. One changes the matrix “in-place” while the other creates a new matrix that is a modified version of the original. This is an important distinction that you should understand for every new Sage command you learn that might change a matrix or vector. Consider the first row operation, which swaps two rows. There are two matrix methods to do this, a “with” version that will create a new, changed matrix, which you will likely want to save, and a plain version that will change the matrix it operates on “in-place.”. The copy() function, which is a general-purpose command, is a way to make a copy of a matrix before you make changes to it. Study the example below carefully, and then read the explanation following. (Remember that counting begins with zero.) A = matrix(QQ,2,3,[1,2,3,4,5,6]) B = A C = copy(A) D = A.with_swapped_rows(0,1) D[0,0] = -1 A.swap_rows(0,1) A[1,2] = -6 A [ 4 5 6] [ 1 2 -6] B [ 4 5 6] [ 1 2 -6] Chapter SLE 10 C [1 2 3] [4 5 6] D [-1 5 6] [ 1 2 3] Here is how each of these four matrices comes into existence. 1. A is our original matrix, created from a list of entries. 2. B is not a new copy of A, it is just a new name for referencing the exact same matrix internally. 3. C is a brand new matrix, stored internally separate from A, but with identical contents. 4. D is also a new matrix, which is created by swapping the rows of A And here is how each matrix is affected by the commands. 1. A is changed twice “in-place”. First, its rows are swapped rows a “plain” matrix method. Then its entry in the lower-right corner is set to -6. 2. B is just another name for A. So whatever changes are made to A will be evident when we ask for the matrix by the name B. And vice-versa. 3. C is a copy of the original A and does not change, since no subsequent com- mands act on C. 4. D is a new copy of A, created by swapping the rows of A. Once created from A, it has a life of its own, as illustrated by the change in its entry in the upper-left corner to -1. An interesting experiment is to rearrange some of the lines above (or add new ones) and predict the result. Just as with the first row operation, swapping rows, Sage supports the other two row operations with natural sounding commands, with both “in-place” versions and new-matrix versions. A = matrix(QQ, 3, 4, [[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]]) B = copy(A) A.rescale_row(1, 1/2) A [ 1 2 3 4] [5/2 3 7/2 4] [ 9 10 11 12] A.add_multiple_of_row(2, 0, 10) A Chapter SLE 11 [ 1 2 3 4] [5/2 3 7/2 4] [ 19 30 41 52] B.with_rescaled_row(1, 1/2) [ 1 2 3 4] [5/2 3 7/2 4] [ 9 10 11 12] C = B.with_added_multiple_of_row(2, 0, 10) C [ 1 2 3 4] [ 5 6 7 8] [19 30 41 52] Notice that the order of the arguments might feel a bit odd, compared to how we write and think about row operations. Also note how the “with” versions leave a trail of matrices for each step while the plain versions just keep changing A. RREF Reduced Row-Echelon Form There has been a lot of information about using Sage with vectors and matrices in this section. But we can now construct the coefficient matrix of a system of equations and the vector of constants. From these pieces we can easily construct the augmented matrix, which we could subject to a series of row operations. Computers are suppose to make routine tasks easy so we can concentrate on bigger ideas. No exception here, Sage can bring a matrix (augmented or not) to reduced row echelon form with no problem. Let us redo Example SAB with Sage. coeff = matrix(QQ, [[-7, -6, -12], [5, 5, 7], [1, 0, 4]]) const = vector(QQ, [-33, 24, 5]) aug = coeff.augment(const) aug.rref() [ 1 0 0 -3] [ 0 1 0 5] [ 0 0 1 2] And the solution x1 = −3, x2 = 5, x3 = 2 is now obvious. Beautiful. You may notice that Sage has some commands with the word “echelon” in them. For now, these should be avoided like the plague, as there are some subtleties in how they work. The matrix method .rref() will be sufficient for our purposes for a long, long time — so stick with it. Section TSS: Types of Solution Sets FDV Free and Dependent Variables Sage has the matrix method .pivot() to quickly and easily identify the pivot columns of the reduced row-echelon form of a matrix. Notice that we do not have to row-reduce the matrix first, we just ask which columns of a matrix A would be Chapter SLE 12 the pivot columns of the matrix B that is row-equivalent to A and in reduced row- echelon form. By Definition IDV, the indices of the pivot columns for an augmented matrix of a system of equations are the indices of the dependent variables. And the remainder are free variables. But be careful, Sage numbers columns starting from zero and mathematicians typically number variables starting from one. Let us reprise Example ISSI. coeff = matrix(QQ, [[ 1, 4, 0, -1, 0, 7, -9], [ 2, 8,-1, 3, 9, -13, 7], [ 0, 0, 2, -3, -4, 12, -8], [-1, -4, 2, 4, 8, -31, 37]]) const = vector(QQ, [3, 9, 1, 4]) aug = coeff.augment(const) dependent = aug.pivots() dependent (0, 2, 3) So, incrementing each column index by 1 gives us the same set D of indices for the dependent variables. To get the free variables, we can use the following code. Study it and then read the explanation following. free = [index for index in range(7) if not index in dependent] free [1, 4, 5, 6] This is a Python programming construction known as a “list comprehension” but in this setting I prefer to call it “set builder notation.” Let us dissect the command in pieces. The brackets ([,]) create a new list. The items in the list will be values of the variable index. range(7) is another list, integers starting at 0 and stopping just before 7. (While perhaps a bit odd, this works very well when we consistently start counting at zero.) So range(7) is the list [0,1,2,3,4,5,6]. Think of these as candidate values for index, which are generated by for index in range(7). Then we test each candidate, and keep it in the new list if it is not in the list dependent. This is entirely analogous to the following mathematics: F = { f | 1 ≤ f ≤ 7, f 6∈ D} where F is free, f is index, and D is dependent, and we make the 0/1 counting adjustments. This ability to construct sets in Sage with notation so closely mirroring the mathematics is a powerful feature worth mastering. We will use it repeatedly. It was a good exercise to use a list comprehension to form the list of columns that are not pivot columns. However, Sage has us covered. free_and_easy = coeff.nonpivots() free_and_easy (1, 4, 5, 6) Can you use this new matrix method to make a simpler version of the consistent() function we designed above? Chapter SLE 13 RCLS Recognizing Consistency of a Linear System Another way of expressing Theorem RCLS is to say a system is consistent if and only if column n + 1 is not a pivot column of B. Sage has the matrix method .pivot() to easily identify the pivot columns of a matrix. Let us use Archetype E as an example. coeff = matrix(QQ, [[ 2, 1, 7, -7], [-3, 4, -5, -6], [ 1, 1, 4, -5]]) const = vector(QQ, [2, 3, 2]) aug = coeff.augment(const) aug.rref() [ 1 0 3 -2 0] [ 0 1 1 -3 0] [ 0 0 0 0 1] aug.pivots() (0, 1, 4) We can look at the reduced row-echelon form of the augmented matrix and see a pivot column in the final column, so we know the system is inconsistent. However, we could just as easily not form the reduced row-echelon form and just look at the list of pivot columns computed by aug.pivots(). Since aug has 5 columns, the final column is numbered 4, which is present in the list of pivot columns, as we expect. One feature of Sage is that we can easily extend its capabilities by defining new commands. Here will create a function that checks if an augmented matrix represents a consistent system of equations. The syntax is just a bit complicated. lambda is the word that indicates we are making a new function, the input is tem- porarily named A (think Augmented), and the name of the function is consistent. Everything following the colon will be evaluated and reported back as the output. consistent = lambda A: not(A.ncols()-1 in A.pivots()) Execute this block above. There will not be any output, but now the consistent function will be defined and available. Now give it a try (after making sure to have run the code above that defines aug). Note that the output of consistent() will be either True or False. consistent(aug) False The consistent() command works by simply checking to see if the last column of A is not in the list of pivots. We can now test many different augmented matrices, such as perhaps changing the vector of constants while keeping the coefficient matrix fixed. Again, make sure you execute the code above that defines coeff and const. consistent(coeff.augment(const)) False Chapter SLE 14 w = vector(QQ, [3,1,2]) consistent(coeff.augment(w)) True u = vector(QQ, [1,3,1]) consistent(coeff.augment(u)) False Why do some vectors of constants lead to a consistent system with this coefficient matrix, while others do not? This is a fundamental question, which we will come to understand in several different ways. SS1 Solving Systems, Part 1 Sage has a built-in command that will solve a linear system of equations, given a coefficient matrix and a vector of constants. We need to learn some more theory before we can entirely understand this command, but we can begin to explore its use. For now, consider these methods experimental and do not let it replace row- reducing augmented matrices. The matrix method A.solve_right(b) will provide information about solutions to the linear system of equations with coefficient matrix A and vector of constants b. The reason for the “right” (and the corresponding command named with “left”) will have to wait for Sage MVP. For now, it is generally correct in this course to use the “right” variant of any Sage linear algebra command that has both “left” and “right” variants. Let us apply the .solve_right() command to a system with no solutions, in particular Archetype E. We have already seen in Sage RCLS that this system is inconsistent. coeff = matrix(QQ, [[ 2, 1, 7, -7], [-3, 4, -5, -6], [ 1, 1, 4, -5]]) const = vector(QQ, [2, 3, 2]) coeff.solve_right(const) Traceback (most recent call last): ... ValueError: matrix equation has no solutions This is our first discussion of Sage error messages, though undoubtedly you have seen several already! First, here we only show the first and last lines of the message since typically it contains a lot of information specific to whichever computer you may be using. but we always begin with the last line as the most important indication of what has happened. Here the “problem” is quite evident: we get an “error” message telling us that the matrix equation has no solutions. We can debate whether or not this is really an error, but that is the design decision taken in Sage — we just need to be aware of it, the .solve_right() is really only valuable when there is a solution. Generally, when deciphering Sage error messages, you want to start at the bot- tom of the “traceback” and read up through the various routines that have been called. Execute the block above and you may see references to matrix methods such as ._solve_right_general() and then .solve_right(). With time and practice, these mysterious messages will become more and more helpful, so spend some time reading them in tandem with locating the real source of any problems you encounter. What does .solve_right() do with a system that does have solutions? Let us take a look at Example ISSI again, as we did in Sage FDV. Chapter SLE 15 coeff = matrix(QQ, [[ 1, 4, 0, -1, 0, 7, -9], [ 2, 8, -1, 3, 9, -13, 7], [ 0, 0, 2, -3, -4, 12, -8], [-1, -4, 2, 4, 8, -31, 37]]) const = vector(QQ, [3, 9, 1, 4]) coeff.solve_right(const) (4, 0, 2, 1, 0, 0, 0) This vector with 7 entries is indeed a solution to the system of equations (check this!). But from Example ISSI we know this system has infinitely many solutions. Why does Sage give us just one solution? Of the infinitely many solutions, why this one? How can a finite computer ever present us with infinitely many solutions? Do we have the time to read through an infinite list of solutions? Is there a “best” solution? This behavior should prompt these questions, and maybe more. In order to totally understand the behavior of the .solve_right() command, we need to understand more of the theory of linear algebra. In good time. So for now, .solve_right() is a curiosity we will fully understand soon — specifically in Sage SS2 and Sage SS3. Section HSE: Homogeneous Systems of Equations SHS Solving Homogeneous Systems We can explore homogeneous systems easily in Sage with commands we have already learned. Notably, the zero_vector() constructor will quickly create the necessary vector of constants (Sage V). You could try defining a function to accept a matrix as input, augment the matrix with a zero vector having the right number of entries, and then return the reduced row-echelon form of the augmented matrix, or maybe you could return the number of free variables. (See Sage RCLS for a refresher on how to do this). It is also interesting to see how .solve_right() behaves on homogeneous systems, since in particular we know it will never give us an error. (Why not? Hint: Theorem HSC.) NS Null Space Sage will compute a null space for us. Which is rather remarkable, as it is an infinite set! Again, this is a powerful command, and there is lots of associated theory, so we will not understand everything about it right away, and it also has a radically different name in Sage. But we will find it useful immediately. Let us reprise Example NSEAI. The relevant command to build the null space of a matrix is .right_kernel(), where again, we will rely exclusively on the “right” version. Also, to match our work in the text, and make the results more recognizable, we will consistently us the keyword option basis=’pivot’, which we will be able to explain once we have more theory (Sage SSNS, Sage SUTH0). Note too, that this is a place where it is critical that matrices are defined to use the rationals as their number system (QQ). I = matrix(QQ, [[ 1, 4, 0, -1, 0, 7, -9], [ 2, 8, -1, 3, 9, -13, 7], [ 0, 0, 2, -3, -4, 12, -8], [-1, -4, 2, 4, 8, -31, 37]]) nsp = I.right_kernel(basis=’pivot’) nsp Chapter SLE 16 Vector space of degree 7 and dimension 4 over Rational Field User basis matrix: [-4 1 0 0 0 0 0] [-2 0 -1 -2 1 0 0] [-1 0 3 6 0 1 0] [ 3 0 -5 -6 0 0 1] As we said, nsp contains a lot of unfamiliar information. Ignore most of it for now. But as a set, we can test membership in nsp. x = vector(QQ, [3, 0, -5, -6, 0, 0, 1]) x in nsp True y = vector(QQ, [-4, 1, -3, -2, 1, 1, 1]) y in nsp True z = vector(QQ, [1, 0, 0, 0, 0, 0, 2]) z in nsp False We √ did a bad thing above, as Sage likes to use I for the imaginary number i = −1 and we just clobbered that. We will not do it again. See below how to fix this. nsp is an infinite set. Since we know the null space is defined as solution to a system of equations, and the work above shows it has at least two elements, we are not surprised to discover that the set is infinite (Theorem PSSLS). nsp.is_finite() False If we want an element of the null space to experiment with, we can get a “ran- dom” element easily. Evaluate the following compute cell repeatedly to get a feel for the variety of the different output. You will see a different result each time, and the result supplied in your downloaded worksheet is very unlikely to be a result you will ever see again. The bit of text, # random, is technically a “comment”, but we are using it as a signal to our automatic testing of the Sage examples that this example should be skipped. You do not need to use this device in your own work, though you may use the comment syntax if you wish. z = nsp.random_element() z # random (21/5, 1, -102/5, -204/5, -3/5, -7, 0) z in nsp Chapter SLE 17 True Sometimes, just sometimes, the null space is finite, and we can list its elements. This is from Example CNS2. C = matrix(QQ, [[-4, 6, 1], [-1, 4, 1], [ 5, 6, 7], [ 4, 7, 1]]) Cnsp = C.right_kernel(basis=’pivot’) Cnsp.is_finite() True Cnsp.list() [(0, 0, 0)] Notice that we get back a list (which mathematically is really a set), and it has one element, the three-entry zero vector. There is more to learn about exploring the null space with Sage’s .right_kernel() so we will see more of this matrix method. In the meantime, if you are done experi- menting with√the matrix I we can restore the variable I back to being the imaginary number i = −1 with the Sage restore() command. restore() I^2 -1 SH Sage Help There are many ways to learn about, or remind yourself of, how various Sage com- mands behave. Now that we have learned a few, it is a good time to show you the most direct methods of obtaining help. These work throughout Sage, so can be useful if you want to apply Sage to other areas of mathematics. The first hurdle is to learn how to make a mathematical object in Sage. We know now how to make matrices and vectors (and null spaces). This is enough to help us explore relevant commands in Sage for linear algebra. First, define a very simple matrix A, with maybe one with one row and two columns. The number system you choose will have some effect on the results, so use QQ for now. In the notebook, enter A. (assuming you called your matrix A, and be sure to include the period). Now hit the “tab” key and you will get a long list of all the possible methods you can apply to A using the dot notation. You can click directly on one of these commands (the word, not the blue high- light) to enter it into the cell. Now instead of adding parentheses to the command, place a single question mark (?) on the end and hit the tab key again. You should get some nicely formatted documentation, along with example uses. (Try A.rref? below for a good example of this.) You can replace the single question mark by two question marks, and as Sage is an open source program you can see the actual computer instructions for the method, which at first includes all the documentation again. Note that now the documentation is enclosed in a pair of triple quotation marks (""",""") as part of the source code, and is not specially formatted. These methods of learning about Sage are generally referred to as “tab-completion” and we will use this term going forward. To learn about the use of Sage in other areas of mathematics, you just need to find out how to create the relevant objects Chapter SLE 18 via a “constructor” function, such as matrix() and vector() for linear algebra. Sage has a comprehensive Reference Manual and there is a Linear Algebra Quick Reference sheet. These should be easily located online via sagemath.org or with an internet search leading with the terms “sage math” (use “math” to avoid confusion with other websites for things named “Sage”). Section NM: Nonsingular Matrices NM Nonsingular Matrix Being nonsingular is an important matrix property, and in such cases Sage contains commands that quickly and easily determine if the mathematical object does, or does not, have the property. The names of these types of methods universally begin with .is_, and these might be referred to as “predicates” or “queries.”. In the Sage notebook, define a simple matrix A, and then in a cell type A.is_, followed by pressing the tab key rather than evaluating the cell. You will get a list of numerous properties that you can investigate for the matrix A. (This will not work as advertised with the Sage cell server.) The other convention is to name these properties in a positive way, so the relevant command for nonsingular matrices is .is_singular(). We will redo Example S and Example NM. Note the use of not in the last compute cell. A = matrix(QQ, [[1, -1, 2], [2, 1, 1], [1, 1, 0]]) A.is_singular() True B = matrix(QQ, [[-7, -6, -12], [ 5, 5, 7], [ 1, 0, 4]]) B.is_singular() False not(B.is_singular()) True IM Identity Matrix It is straightforward to create an identity matrix in Sage. Just specify the number system and the number of rows (which will equal the number of columns, so you do not specify that since it would be redundant). The number system can be left out, but the result will have entries from the integers (ZZ), which in this course is unlikely to be what you really want. id5 = identity_matrix(QQ, 5) id5 [1 0 0 0 0] [0 1 0 0 0] [0 0 1 0 0] [0 0 0 1 0] Chapter SLE 19 [0 0 0 0 1] id4 = identity_matrix(4) id4.base_ring() Integer Ring Notice that we do not use the now-familiar dot notation to create an identity matrix. What would we use the dot notation on anyway? For these reasons we call the identity_matrix() function a constructor, since it builds something from scratch, in this case a very particular type of matrix. We mentioned above that an identity matrix is in reduced row-echelon form. What happens if we try to row-reduce a matrix that is already in reduced row- echelon form? By the uniqueness of the result, there should be no change. The following code illustrates this. Notice that = is used to assign an object to a new name, while == is used to test equality of two objects. I frequently make the mistake of forgetting the second equal sign when I mean to test equality. id50 = identity_matrix(QQ, 50) id50 == id50.rref() True NME1 Nonsingular Matrix Equivalences, Round 1 Sage will create random matrices and vectors, sometimes with various properties. These can be very useful for quick experiments, and they are also useful for illus- trating that theorems hold for any object satisfying the hypotheses of the theorem. But this will never replace a proof. We will illustrate Theorem NME1 using Sage. We will use a variant of the random_matrix() constructor that uses the algorithm=’unimodular’ keyword. We will have to wait for Chapter D before we can give a full explanation, but for now, understand that this command will always create a square matrix that is nonsingular. Also realize that there are square nonsingular matrices which will never be the output of this command. In other words, this command creates elements of just a subset of all possible nonsingular matrices. So we are using random matrices below to illustrate properties predicted by Theorem NME1. Execute the first command to create a random nonsingular matrix, and notice that we only have to mark the output of A as random for our automated testing process. After a few runs, notice that you can also edit the value of n to create matrices of different sizes. With a matrix A defined, run the next three cells, which by Theorem NME1 each always produce True as their output, no matter what value A has, so long as A is nonsingular. Read the code and try to determine exactly how they correspond to the parts of the theorem (some commentary along these lines follows). n = 6 A = random_matrix(QQ, n, algorithm=’unimodular’) A # random [ 1 -4 8 14 8 55] [ 4 -15 29 50 30 203] [ -4 17 -34 -59 -35 -235] [ -1 3 -8 -16 -5 -48] [ -5 16 -33 -66 -16 -195] Chapter SLE 20 [ 1 -2 2 7 -2 10] A.rref() == identity_matrix(QQ, n) True nsp = A.right_kernel(basis=’pivot’) nsp.list() == [zero_vector(QQ, n)] True b = random_vector(QQ, n) aug = A.augment(b) aug.pivots() == tuple(range(n)) True The only portion of these commands that may be unfamilar is the last one. The command range(n) is incredibly useful, as it will create a list of the integers from 0 up to, but not including, n. (We saw this command briefly in Sage FDV.) So, for example, range(3) == [0,1,2] is True. Pivots are returned as a “tuple” which is very much like a list, except we cannot change the contents. We can see the difference by the way the tuple prints with parentheses ((,)) rather than brackets ([,]). We can convert a list to a tuple with the tuple() command, in order to make the comparison succeed. How do we tell if the reduced row-echelon form of the augmented matrix of a system of equations represents a system with a unique solution? First, the system must be consistent, which by Theorem RCLS means the last column is not a pivot column. Then with a consistent system we need to insure there are no free variables. This happens if and only if the remaining columns are all pivot columns, according to Theorem FVCS, thus the test used in the last compute cell. Chapter V Vectors Section VO: Vector Operations VSCV Vector Spaces of Column Vectors It is possible to construct vector spaces several ways in Sage. For now, we will show you two basic ways. Remember that while our theory is all developed over the complex numbers, C, it is better to initially illustrate these ideas in Sage using the rationals, QQ. To create a vector space, we use the VectorSpace() constructor, which requires the name of the number system for the entries and the number of entries in each vector. We can display some information about the vector space, and with tab- completion you can see what functions are available. We will not do too much with these methods immediately, but instead learn about them as we progress through the theory. V = VectorSpace(QQ, 8) V Vector space of dimension 8 over Rational Field Notice that the word “dimension” is used to refer to the number of entries in a vector contained in the vector space, whereas we have used the word “degree” before. Try pressing the Tab key while in the next cell to see the range of methods you can use on a vector space. V. We can easily create “random” elements of any vector space, much as we did earlier for the kernel of a matrix. Try executing the next compute cell several times. w = V.random_element() w # random (2, -1/9, 0, 2, 2/3, 0, -1/3, 1) Vector spaces are a fundamental objects in Sage and in mathematics, and Sage has a nice compact way to create them, mimicking the notation we use when working on paper. U = CC^5 U 21 Chapter V 22 Vector space of dimension 5 over Complex Field with 53 bits of precision W = QQ^3 W Vector space of dimension 3 over Rational Field Sage can determine if two vector spaces are the same. Notice that we use two equals sign to test equality, since we use a single equals sign to make assignments. X = VectorSpace(QQ, 3) W = QQ^3 X == W True VO Vector Operations Sage can easily perform the two basic operations with vectors, vector addition, +, and scalar vector multiplication, *. Notice that Sage is not confused by an ambiguity due to multiple meanings for the symbols + and * — for example, Sage knows that 3 + 12 is different than the vector additions below. x = vector(QQ, [1, 2, 3]) y = vector(QQ, [10, 20, 30]) 5*x (5, 10, 15) x + y (11, 22, 33) 3*x + 4*y (43, 86, 129) -y (-10, -20, -30) w = (-4/3)*x - (1/10)*y w (-7/3, -14/3, -7) ANC A Note on Coercion Study the following sequence of commands, while cognizant of the failure to specify a number system for x. Chapter V 23 x = vector([1, 2, 3]) u = 3*x u (3, 6, 9) v = (1/3)*x v (1/3, 2/3, 1) y = vector(QQ, [4, 5, 6]) w = 8*y w (32, 40, 48) z = x + y z (5, 7, 9) None of this should be too much of a surprise, and the results should be what we would have expected. Though for x we never specified if 1, 2, 3 are integers, rationals, reals, complexes, or . . . ? Let us dig a little deeper and examine the parents of the five vectors involved. x.parent() Ambient free module of rank 3 over the principal ideal domain Integer Ring u.parent() Ambient free module of rank 3 over the principal ideal domain Integer Ring v.parent() Vector space of dimension 3 over Rational Field y.parent() Vector space of dimension 3 over Rational Field Chapter V 24 w.parent() Vector space of dimension 3 over Rational Field z.parent() Vector space of dimension 3 over Rational Field So x and u belong to something called an “ambient free module,” whatever that is. What is important here is that the parent of x uses the integers as its number system. How about u, v, y, w, z? All but the first has a parent that uses the rationals for its number system. Three of the final four vectors are examples of a process that Sage calls “coer- cion.” Mathematical elements get converted to a new parent, as necessary, when the conversion is totally unambiguous. In the examples above: • u is the result of scalar multiplication by an integer, so the computation and result can all be accommodated within the integers as the number system. • v involves scalar multiplication by a scalar that is not an integer, and which could be construed as a rational number. So the result needs to have a parent whose number system is the rationals. • y is created explicitly as a vector whose entries are rational numbers. • Even though w is created only with products of integers, the fact that y has entries considered as rational numbers, so too does the result. • The creation of z is the result of adding a vector of integers to a vector of rationals. This is the best example of coercion — Sage promotes x to a vector of rationals and therefore returns a result that is a vector of rationals. Notice that there is no ambiguity and no argument about how to promote x, and the same would be true for any vector full of integers. The coercion above is automatic, but we can also usually force it to happen without employing an operation. t = vector([10, 20, 30]) t.parent() Ambient free module of rank 3 over the principal ideal domain Integer Ring V = QQ^3 t_rational = V(t) t_rational (10, 20, 30) t_rational.parent() Vector space of dimension 3 over Rational Field Chapter V 25 W = CC^3 t_complex = W(t) t_complex (10.0000000000000, 20.0000000000000, 30.0000000000000) t_complex.parent() Vector space of dimension 3 over Complex Field with 53 bits of precision So the syntax is to use the name of the parent like a function and coerce the element into the new parent. This can fail if there is no natural way to make the conversion. u = vector(CC, [5*I, 4-I]) u (5.00000000000000*I, 4.00000000000000 - 1.00000000000000*I) V = QQ^2 V(u) Traceback (most recent call last): ... TypeError: Unable to coerce 5.00000000000000*I (<type ’sage.rings.complex_number.ComplexNumber’>) to Rational Coercion is one of the more mysterious aspects of Sage, and the above discussion may not be very clear the first time though. But if you get an error (like the one above) talking about coercion, you know to come back here and have another read through. For now, be sure to create all your vectors and matrices over QQ and you should not have any difficulties. Section LC: Linear Combinations LC Linear Combinations We can redo Example TLC with Sage. First we build the relevant vectors and then do the computation. u1 = vector(QQ, [ 2, 4, -3, 1, 2, 9]) u2 = vector(QQ, [ 6, 3, 0, -2, 1, 4]) u3 = vector(QQ, [-5, 2, 1, 1, -3, 0]) u4 = vector(QQ, [ 3, 2, -5, 7, 1, 3]) 1*u1 + (-4)*u2 + 2*u3 +(-1)*u4 (-35, -6, 4, 4, -9, -10) With a linear combination combining many vectors, we sometimes will use more compact ways of forming a linear combination. So we will redo the second linear combination of u1 , u2 , u3 , u4 using a list comprehension and the sum() function. Chapter V 26 vectors = [u1, u2, u3, u4] scalars = [3, 0, 5, -1] multiples = [scalars[i]*vectors[i] for i in range(4)] multiples [(6, 12, -9, 3, 6, 27), (0, 0, 0, 0, 0, 0), (-25, 10, 5, 5, -15, 0), (-3, -2, 5, -7, -1, -3)] We have constructed two lists and used a list comprehension to just form the scalar multiple of each vector as part of the list multiples. Now we use the sum() function to add them all together. sum(multiples) (-22, 20, 1, 1, -10, 24) We can improve on this in two ways. First, we can determine the number of elements in any list with the len() function. So we do not have to count up that we have 4 vectors (not that it is very hard to count!). Second, we can combine this all into one line, once we have defined the list of vectors and the list of scalars. sum([scalars[i]*vectors[i] for i in range(len(vectors))]) (-22, 20, 1, 1, -10, 24) The corresponding expression in mathematical notation, after a change of names and with counting starting from 1, would roughly be: 4 X ai ui i=1 Using sum() and a list comprehension might be overkill in this example, but we will find it very useful in just a minute. SLC Solutions and Linear Combinations We can easily illustrate Theorem SLSLC with Sage. We will use Archetype F as an example. coeff = matrix(QQ, [[33, -16, 10,-2], [99, -47, 27,-7], [78, -36, 17,-6], [-9, 2, 3, 4]]) const = vector(QQ, [-27, -77, -52, 5]) A solution to this system is x1 = 1, x2 = 2, x3 = −2, x4 = 4. So we will use these four values as scalars in a linear combination of the columns of the coefficient matrix. However, we do not have to type in the columns individually, we can have Sage extract them all for us into a list with the matrix method .columns(). cols = coeff.columns() cols Chapter V 27 [(33, 99, 78, -9), (-16, -47, -36, 2), (10, 27, 17, 3), (-2, -7, -6, 4)] With our scalars also in a list, we can compute the linear combination of the columns, like we did in Sage LC. soln = [1, 2, -2, 4] sum([soln[i]*cols[i] for i in range(len(cols))]) (-27, -77, -52, 5) So we see that the solution gives us scalars that yield the vector of constants as a linear combination of the columns of the coefficient matrix. Exactly as predicted by Theorem SLSLC. We can duplicate this observation with just one line: const == sum([soln[i]*cols[i] for i in range(len(cols))]) True In a similar fashion we can test other potential solutions. With theory we will develop later, we will be able to determine that Archetype F has only one solution. Since Theorem SLSLC is an equivalence (Technique E), any other choice for the scalars should not create the vector of constants as a linear combination. alt_soln = [-3, 2, 4, 1] const == sum([alt_soln[i]*cols[i] for i in range(len(cols))]) False Now would be a good time to find another system of equations, perhaps one with infinitely many solutions, and practice the techniques above. SS2 Solving Systems, Part 2 We can now resolve a bit of the mystery around Sage’s .solve_right() method. Recall from Sage SS1 that if a linear system has solutions, Sage only provides one solution, even in the case when there are infinitely many solutions. In our previous discussion, we used the system from Example ISSI. coeff = matrix(QQ, [[ 1, 4, 0, -1, 0, 7, -9], [ 2, 8, -1, 3, 9, -13, 7], [ 0, 0, 2, -3, -4, 12, -8], [-1, -4, 2, 4, 8, -31, 37]]) const = vector(QQ, [3, 9, 1, 4]) coeff.solve_right(const) (4, 0, 2, 1, 0, 0, 0) The vector c described in the statement of Theorem VFSLS is precisely the solution returned from Sage’s .solve_right() method. This is the solution where we choose the αi , 1 ≤ i ≤ n−r to all be zero, in other words, each free variable is set to zero (how convenient!). Free variables correspond to columns of the row-reduced augmented matrix that are not pivot columns. So we can profitably employ the .nonpivots() matrix method. Let us put this all together. Chapter V 28 aug = coeff.augment(const) reduced = aug.rref() reduced [ 1 4 0 0 2 1 -3 4] [ 0 0 1 0 1 -3 5 2] [ 0 0 0 1 2 -6 6 1] [ 0 0 0 0 0 0 0 0] aug.nonpivots() (1, 4, 5, 6, 7) Since the eighth column (numbered 7) of the reduced row-echelon form is not a pivot column, we know by Theorem RCLS that the system is consistent. We can use the indices of the remaining non-pivot columns to place zeros into the vector c in those locations. The remaining entries of c are the entries of the reduced row-echelon form in the last column, inserted in that order. Boom! So we have three ways to get to the same solution: (a) row-reduce the aug- mented matrix and set the free variables all to zero, (b) row-reduce the augmented matrix and use the formula from Theorem VFSLS to construct c, and (c) use Sage’s .solve_right() method. One mystery left to resolve. How can we get Sage to give us infinitely many solutions in the case of systems with an infinite solution set? This is best handled in the next section, Section SS, specifically in Sage SS3. PSHS Particular Solutions, Homogeneous Solutions Again, Sage is useful for illustrating a theorem, in this case Theorem PSPHS. We will illustrate both “directions” of this equivalence with the system from Example ISSI. coeff = matrix(QQ,[[ 1, 4, 0, -1, 0, 7, -9], [ 2, 8, -1, 3, 9, -13, 7], [ 0, 0, 2, -3, -4, 12, -8], [-1, -4, 2, 4, 8, -31, 37]]) n = coeff.ncols() const = vector(QQ, [3, 9, 1, 4]) First we will build solutions to this system. Theorem PSPHS says we need a particular solution, i.e. one solution to the system, w. We can get this from Sage’s .solve_right() matrix method. Then for any vector z from the null space of the coefficient matrix, the new vector y = w + z should be a solution. We walk through this construction in the next few cells, where we have employed a specific element of the null space, z, along with a check that it is really in the null space. w = coeff.solve_right(const) nsp = coeff.right_kernel(basis=’pivot’) z = vector(QQ, [42, 0, 84, 28, -50, -47, -35]) z in nsp True Chapter V 29 y = w + z y (46, 0, 86, 29, -50, -47, -35) const == sum([y[i]*coeff.column(i) for i in range(n)]) True You can create solutions repeatedly via the creation of random elements of the null space. Be sure you have executed the cells above, so that coeff, n, const, nsp and w are all defined. Try executing the cell below repeatedly to test infinitely many solutions to the system. You can use the subsequent compute cell to peek at any of the solutions you create. z = nsp.random_element() y = w + z const == sum([y[i]*coeff.column(i) for i in range(n)]) True y # random (-11/2, 0, 45/2, 34, 0, 7/2, -2) For the other direction, we present (and verify) two solutions to the linear system of equations. The condition that y = w + z can be rewritten as y − w = z, where z is in the null space of the coefficient matrix. which of our two solutions is the “particular” solution and which is “some other” solution? It does not matter, it is all sematics at this point. What is important is that their difference is an element of the null space (in either order). So we define the solutions, along with checks that they are really solutions, then examine their difference. soln1 = vector(QQ,[4, 0, -96, 29, 46, 76, 56]) const == sum([soln1[i]*coeff.column(i) for i in range(n)]) True soln2 = vector(QQ,[-108, -84, 86, 589, 240, 283, 105]) const == sum([soln2[i]*coeff.column(i) for i in range(n)]) True (soln1 - soln2) in nsp True Chapter V 30 Section SS: Spanning Sets SS Spanning Sets A strength of Sage is the ability to create infinite sets, such as the span of a set of vectors, from finite descriptions. In other words, we can take a finite set with just a handful of vectors and Sage will create the set that is the span of these vectors, which is an infinite set. Here we will show you how to do this, and show how you can use the results. The key command is the vector space method .span(). V = QQ^4 v1 = vector(QQ, [1,1,2,-1]) v2 = vector(QQ, [2,3,5,-4]) W = V.span([v1, v2]) W Vector space of degree 4 and dimension 2 over Rational Field Basis matrix: [ 1 0 1 1] [ 0 1 1 -2] x = 2*v1 + (-3)*v2 x (-4, -7, -11, 10) x in W True y = vector(QQ, [3, -1, 2, 2]) y in W False u = vector(QQ, [3, -1, 2, 5]) u in W True W <= V True Most of the above should be fairly self-explanatory, but a few comments are in order. The span, W, is created in the first compute cell with the .span() method, which accepts a list of vectors and needs to be employed as a method of a vector space. The information about W printed when we just input the span itself may be somewhat confusing, and as before, we need to learn some more theory to totally understand it all. For now you can see the number system (Rational Field) and the number of entries in each vector (degree 4). The dimension may be more complicated than you first suspect. Chapter V 31 Sets are all about membership, and we see that we can easily check membership in a span. We know the vector x will be in the span W since we built it as a linear combination of v1 and v2. The vectors y and u are a bit more mysterious, but Sage can answer the membership question easily for both. The last compute cell is something new. The symbol <= is meant here to be the “subset of” relation, i.e. a slight abuse of the mathematical symbol ⊆, and then we are not surprised to learn that W really is a subset of V. It is important to realize that the span is a construction that begins with a finite set, yet creates an infinite set. With a loop (the for command) and the .random_element() vector space method we can create many, but not all, of the elements of a span. In the examples below, you can edit the total number of random vectors produced, and you may need to click through to another file of output if you ask for more than about 100 vectors. Each example is designed to illustrate some aspect of the behavior of the span command and to provoke some questions. So put on your mathematician’s hat, evaluate the compute cells to create some sample elements, and then study the output carefully looking for patterns and maybe even conjecture some explanations for the behavior. The puzzle gets just a bit harder for each new example. (We use the Sequence() command to get nicely-formatted line-by-line output of the list, and notice that we are only providing a portion of the output here. You will want to evalaute the computation of vecs and then evaluate the next cell in the Sage notebook for maximum effect.) V = QQ^4 W = V.span([ vector(QQ, [0, 1, 0, 1]), vector(QQ, [1, 0, 1, 0]) ]) vecs = [(i, W.random_element()) for i in range(100)] Sequence(vecs, cr=True) # random [ (0, (1/5, 16, 1/5, 16)), (1, (-3, 0, -3, 0)), (2, (1/11, 0, 1/11, 0)), ... (97, (-2, -1/2, -2, -1/2)), (98, (1/13, -3, 1/13, -3)), (99, (0, 1, 0, 1)) ] V = QQ^4 W = V.span([ vector(QQ, [0, 1, 1, 0]), vector(QQ, [0, 0, 1, 1]) ]) vecs = [(i, W.random_element()) for i in range(100)] Sequence(vecs, cr=True) # random [ (0, (0, 1/9, 2, 17/9)), (1, (0, -1/8, 3/2, 13/8)), (2, (0, 1/2, -1, -3/2)), ... (97, (0, 4/7, 24, 164/7)), (98, (0, -5/2, 0, 5/2)), (99, (0, 13/2, 1, -11/2)) Chapter V 32 ] V = QQ^4 W = V.span([ vector(QQ, [2, 1, 2, 1]), vector(QQ, [4, 2, 4, 2]) ]) vecs = [(i, W.random_element()) for i in range(100)] Sequence(vecs, cr=True) # random [ (0, (1, 1/2, 1, 1/2)), (1, (-1, -1/2, -1, -1/2)), (2, (-1/7, -1/14, -1/7, -1/14)), ... (97, (1/3, 1/6, 1/3, 1/6)), (98, (0, 0, 0, 0)), (99, (-11, -11/2, -11, -11/2)) ] V = QQ^4 W = V.span([ vector(QQ, [1, 0, 0, 0]), vector(QQ, [0, 1 ,0, 0]), vector(QQ, [0, 0, 1, 0]), vector(QQ, [0, 0, 0, 1]) ]) vecs = [(i, W.random_element()) for i in range(100)] Sequence(vecs, cr=True) # random [ (0, (-7/4, -2, -1, 63)), (1, (6, -2, -1/2, -28)), (2, (5, -2, -2, -1)), ... (97, (1, -1/2, -2, -1)), (98, (-1/12, -4, 2, 1)), (99, (2/3, 0, -4, -1)) ] V = QQ^4 W = V.span([ vector(QQ, [1, 2, 3, 4]), vector(QQ, [0,-1, -1, 0]), vector(QQ, [1, 1, 2, 4]) ]) vecs = [(i, W.random_element()) for i in range(100)] Sequence(vecs, cr=True) # random [ (0, (-1, 3, 2, -4)), (1, (-1/27, -1, -28/27, -4/27)), (2, (-7/11, -1, -18/11, -28/11)), ... (97, (1/3, -1/7, 4/21, 4/3)), Chapter V 33 (98, (-1, -14, -15, -4)), (99, (0, -2/7, -2/7, 0)) ] CSS Consistent Systems and Spans With the notion of a span, we can expand our techniques for checking the consistency of a linear system. Theorem SLSLC tells us a system is consistent if and only if the vector of constants is a linear combination of the columns of the coefficient matrix. This is because Theorem SLSLC says that any solution to the system will provide a linear combination of the columns of the coefficient that equals the vector of constants. So consistency of a system is equivalent to the membership of the vector of constants in the span of the columns of the coefficient matrix. Read that last sentence again carefully. We will see this idea again, but more formally, in Theorem CSCS. We will reprise Sage SLC, which is based on Archetype F. We again make use of the matrix method .columns() to get all of the columns into a list at once. coeff = matrix(QQ, [[33, -16, 10, -2], [99, -47, 27, -7], [78, -36, 17, -6], [-9, 2, 3, 4]]) column_list = coeff.columns() column_list [(33, 99, 78, -9), (-16, -47, -36, 2), (10, 27, 17, 3), (-2, -7, -6, 4)] span = (QQ^4).span(column_list) const = vector(QQ, [-27, -77, -52, 5]) const in span True You could try to find an example of a vector of constants which would create an inconsistent system with this coefficient matrix. But there is no such thing. Here is why — the null space of coeff is trivial, just the zero vector. nsp = coeff.right_kernel(basis=’pivot’) nsp.list() [(0, 0, 0, 0)] The system is consistent, as we have shown, so we can apply Theorem PSPHS. We can read Theorem PSPHS as saying any two different solutions of the system will differ by an element of the null space, and the only possibility for this null space vector is just the zero vector. In other words, any two solutions cannot be different. SSNS Spanning Sets for Null Spaces We have seen that we can create a null space in Sage with the .right_kernel() method for matrices. We use the optional argument basis=’pivot’, so we get exactly the spanning set {z1 , z2 , z3 , . . . , zn−r } described in Theorem SSNS. This optional argument will overide Sage’s default spanning set, whose purpose we will explain fully in Sage SUTH0. Here is Example SSNS again, along with a few extras we will comment on afterwards. Chapter V 34 A = matrix(QQ, [[ 1, 3, 3, -1, -5], [ 2, 5, 7, 1, 1], [ 1, 1, 5, 1, 5], [-1, -4, -2, 0, 4]]) nsp = A.right_kernel(basis=’pivot’) N = nsp.basis() N [ (-6, 1, 1, 0, 0), (-4, 2, 0, -3, 1) ] z1 = N[0] z2 = N[1] z = 4*z1 +(-3)*z2 z (-12, -2, 4, 9, -3) z in nsp True sum([z[i]*A.column(i) for i in range(A.ncols())]) (0, 0, 0, 0) We built the null space as nsp, and then asked for its .basis(). For now, a “basis” will give us a spanning set, and with more theory we will understand it better. This is a set of vectors that form a spanning set for the null space, and with the basis=’pivot’ argument we have asked for the spanning set in the format described in Theorem SSNS. The spanning set N is a list of vectors, which we have extracted and named as z1 and z2. The linear combination of these two vectors, named z, will also be in the null space since N is a spanning set for nsp. As verification, we have have used the five entries of z in a linear combination of the columns of A, yielding the zero vector (with four entries) as we expect. We can also just ask Sage if z is in nsp: z in nsp True Now is an appropriate time to comment on how Sage displays a null space when we just ask about it alone. Just as with a span, the number system and the number of entries are easy to see. Again, dimension should wait for a bit. But you will notice now that the Basis matrix has been replaced by User basis matrix. This is a consequence of our request for something other than the default basis (the ’pivot’ basis). As part of its standard information about a null space, or a span, Sage spits out the basis matrix. This is a matrix, whose rows are vectors in a spanning set. This matrix can be requested by itself, using the method .basis_matrix(). It is important to notice that this is very different than the output of .basis() which is a list of vectors. The two objects print very similarly, but even this is different — compare the organization of the brackets and parentheses. Finally the last command should print true for any span or null space Sage creates. If you rerun the commands below, be sure the null space nsp is defined from the code just above. Chapter V 35 nsp Vector space of degree 5 and dimension 2 over Rational Field User basis matrix: [-6 1 1 0 0] [-4 2 0 -3 1] nsp.basis_matrix() [-6 1 1 0 0] [-4 2 0 -3 1] nsp.basis() [ (-6, 1, 1, 0, 0), (-4, 2, 0, -3, 1) ] nsp.basis() == nsp.basis_matrix().rows() True SS3 Solving Systems, Part 3 We have deferred several mysteries and postponed some explanations of the terms Sage uses to describe certain objects. This is because Sage is a comprehensive tool that you can use throughout your mathematical career, and Sage assumes you already know a lot of linear algebra. Do not worry, you already know some linear algebra and will very soon will know a whole lot more. And we know enough now that we can now solve one mystery. How can we create all of the solutions to a linear system of equations when the solution set is infinite? Theorem VFSLS described elements of a solution set as a single vector c, plus linear combinations of vectors u1 , u2 , u3 , . . . , un−r . The vectors uj are identical to the vectors zj in the description of the spanning set of a null space in Theo- rem SSNS, suitably adjusted from the setting of a general system and the relevant homogeneous system for a null space. We can get the single vector c from the .solve_right() method of a coefficient matrix, once supplied with the vector of constants. The vectors {z1 , z2 , z3 , . . . , zn−r } come from Sage in the basis returned by the .right_kernel(basis=’pivot’) method applied to the coefficient matrix. Theorem PSPHS amplifies and generalizes this discussion, making it clear that the choice of the particular particular solution c (Sage’s choice, and your author’s choice) is just a matter of convenience. In our previous discussion, we used the system from Example ISSI. We will use it one more time. coeff = matrix(QQ, [[ 1, 4, 0, -1, 0, 7, -9], [ 2, 8, -1, 3, 9, -13, 7], [ 0, 0, 2, -3, -4, 12, -8], [-1, -4, 2, 4, 8, -31, 37]]) n = coeff.ncols() const = vector(QQ, [3, 9, 1, 4]) c = coeff.solve_right(const) Chapter V 36 c (4, 0, 2, 1, 0, 0, 0) N = coeff.right_kernel(basis=’pivot’).basis() z1 = N[0] z2 = N[1] z3 = N[2] z4 = N[3] soln = c + 2*z1 - 3*z2 - 5*z3 + 8*z4 soln (31, 2, -50, -71, -3, -5, 8) check = sum([soln[i]*coeff.column(i) for i in range(n)]) check (3, 9, 1, 4) check == const True So we can describe the infinitely many solutions with just five items: the specific solution, c, and the four vectors of the spanning set. Boom! As an exercise in understanding Theorem PSPHS, you might begin with soln and add a linear combination of the spanning set for the null space to create another solution (which you can check). Section LI: Linear Independence LI Linear Independence We can use the Sage tools we already have, along with Theorem LIVRN, to de- termine if sets are linearly independent. There is just one hitch — Sage has a preference for placing vectors into matrices as rows, rather than as columns. When printing vectors on the screen, writing them across the screen makes good sense, and there are more mathematical reasons for this choice. But we have chosen to present introductory linear algebra with an emphasis on the columns of matrices — again, for good mathematical reasons. Fortunately, Sage allows us to build matrices from columns as well as rows. Let us redo Example LDHS, the determination that a set of three vectors is linearly dependent. We will enter the vectors, construct a matrix with the vectors as columns via the column_matrix() constructor, and analyze. Here we go. v1 = vector(QQ, [2, -1, 3, 4, 2]) v2 = vector(QQ, [6, 2, -1, 3, 4]) v3 = vector(QQ, [4, 3, -4, -1, 2]) A = column_matrix([v1, v2, v3]) A Chapter V 37 [ 2 6 4] [-1 2 3] [ 3 -1 -4] [ 4 3 -1] [ 2 4 2] A.ncols() == len(A.pivots()) False Notice that we never explicitly row-reduce A, though this computation must happen behind the scenes when we compute the list of pivot columns. We do not really care where the pivots are (the actual list), but rather we want to know how many there are, thus we ask about the length of the list with the function len(). Once we construct the matrix, the analysis is quick. With n 6= r, Theorem LIVRN tells us the set is linearly dependent. Reprising Example LIS with Sage would be good practice at this point. Here is an empty compute cell to use. While it is extremely important to understand the approach outlined above, Sage has a convenient tool for working with linear independence. .linear_dependence() is a method for vector spaces, which we feed in a list of vectors. The output is again a list of vectors, each one containing the scalars that yield a nontrivial relation of linear dependence on the input vectors. We will give this method a workout in the next section, but for now we are interested in the case of a linearly independent set. In this instance, the method will return nothing (an empty list, really). Not even the all-zero vector is produced, since it is not interesting and definitely is not surprising. Again, we will not say anymore about the output of this method until the next section, and do not let its use replace a good conceptual understanding of this section. We will redo Example LDHS again, you try Example LIS again. If you are playing along, be sure v1, v2, v3 are defined from the code above. L = [v1, v2, v3] V = QQ^5 V.linear_dependence(L) == [] False The only comment to make here is that we need to create the vector space QQ^5 since .linear_dependence() is a method of vector spaces. Example LIS should proceed similarly, though being a linearly independent set, the comparison with the empty list should yield True. NME2 Nonsingular Matrix Equivalences, Round 2 As will be our habit, we can illustrate properties of nonsingular matrices with random square matrices. Review Sage NME1 for a refresher on generating these matrices. Here we will illustrate the fifth condition, half of Theorem NMLIC. n = 6 A = random_matrix(QQ, n, algorithm=’unimodular’) A # random [ 2 7 37 -79 268 44] [ -1 -3 -16 33 -110 -16] [ 1 1 7 -14 44 5] [ 0 -3 -15 34 -118 -23] Chapter V 38 [ 2 6 33 -68 227 34] [ -1 -3 -16 35 -120 -21] V = QQ^n V.linear_dependence(A.columns()) == [] True You should always get an empty list ([]) as the result of the second compute cell, no matter which random (unimodular) matrix you produce prior. Note that the list of vectors created using .columns() is exactly the list we want to feed to .linear_dependence(). LISS Linearly Independent Spanning Sets No new commands here, but instead we can now reveal and explore one of our Sage mysteries. When we create a null space with .right_kernel(basis=’pivot’), or a span with .span(), Sage internally creates a new spanning set for the null space or span. We see the vectors of these spanning sets as the rows of the Basis matrix or User basis matrix when we print one of these sets. But more than a spanning set, these vectors are always a linearly independent set. The advantages of this will have to wait for more theory, but for now, this may explain why the spanning set changes (and sometimes shrinks). We can also demonstrate this behavior, by creating the span of a large set of (linearly dependent) vectors. V = QQ^5 v1 = vector(QQ, [ 5, 2, -11, -24, -17]) v2 = vector(QQ, [ 4, 13, -5, -4, 13]) v3 = vector(QQ, [ 8, 29, -9, -4, 33]) v4 = vector(QQ, [ 1, 1, -2, -4, -2]) v5 = vector(QQ, [-7, -25, 8, 4, -28]) v6 = vector(QQ, [ 0, -3, -1, -4, -7]) v7 = vector(QQ, [-1, 2, 3, 8, 9]) L = [v1, v2, v3, v4, v5, v6, v7] V.linear_dependence(L) == [] False W = V.span(L) W Vector space of degree 5 and dimension 2 over Rational Field Basis matrix: [ 1 0 -7/3 -16/3 -13/3] [ 0 1 1/3 4/3 7/3] S = W.basis() S [ (1, 0, -7/3, -16/3, -13/3), (0, 1, 1/3, 4/3, 7/3) ] Chapter V 39 X = V.span(S) X == W True V.linear_dependence(S) == [] True So the above examines properties of the two vectors, S, that Sage computes for the span W built from the seven vectors in the set L. We see that the span of S is equal to the span of L via the comparison X == W. And the smaller set, S, is linearly independent, as promised. Notice that we do not need to use Sage to know that L is linearly dependent, this is guaranteed by Theorem MVSLD. Section LDS: Linear Dependence and Spans RLD Relations of Linear Dependence Example RSC5 turned on a nontrivial relation of linear dependence (Definition RLDCV) on the set {v1 , v2 , v3 , v4 }. Besides indicating linear independence, the Sage vector space method .linear_dependence() produces relations of linear de- pendence for linearly dependent sets. Here is how we would employ this method in Example RSC5. The optional argument zeros=’right’ will produce results con- sistent with our work here, you can also experiment with zeros=’left’ (which is the default). V = QQ^5 v1 = vector(QQ, [1, 2, -1, 3, 2]) v2 = vector(QQ, [2, 1, 3, 1, 2]) v3 = vector(QQ, [0, -7, 6, -11, -2]) v4 = vector(QQ, [4, 1, 2, 1, 6]) R = [v1, v2, v3, v4] L = V.linear_dependence(R, zeros=’right’) L[0] (-4, 0, -1, 1) -4*v1 + 0*v2 +(-1)*v3 +1*v4 (0, 0, 0, 0, 0) V.span(R) == V.span([v1, v2, v4]) True You can check that the list L has just one element (maybe with len(L)), but realize that any multiple of the vector L[0] is also a relation of linear dependence on R, most of which are nontrivial. Notice that we have verified the final conclusion of Example RSC5 with a comparison of two spans. We will give the .linear_dependence() method a real workout in the next Sage subsection (Sage COV) — this is just a quick introduction. Chapter V 40 COV Casting Out Vectors We will redo Example COV, though somewhat tersely, just producing the justifica- tion for each time we toss a vector (a specific relation of linear dependence), and then verifying that the resulting spans, each with one fewer vector, still produce the original span. We also introduce the .remove() method for lists. Ready? Here we go. V = QQ^4 v1 = vector(QQ, [ 1, 2, 0, -1]) v2 = vector(QQ, [ 4, 8, 0, -4]) v3 = vector(QQ, [ 0, -1, 2, 2]) v4 = vector(QQ, [-1, 3, -3, 4]) v5 = vector(QQ, [ 0, 9, -4, 8]) v6 = vector(QQ, [ 7, -13, 12, -31]) v7 = vector(QQ, [-9, 7, -8, 37]) S = [v1, v2, v3, v4, v5, v6, v7] W = V.span(S) D = V.linear_dependence(S, zeros=’right’) D [ (-4, 1, 0, 0, 0, 0, 0), (-2, 0, -1, -2, 1, 0, 0), (-1, 0, 3, 6, 0, 1, 0), (3, 0, -5, -6, 0, 0, 1) ] D[0] (-4, 1, 0, 0, 0, 0, 0) S.remove(v2) W == V.span(S) True D[1] (-2, 0, -1, -2, 1, 0, 0) S.remove(v5) W == V.span(S) True D[2] Chapter V 41 (-1, 0, 3, 6, 0, 1, 0) S.remove(v6) W == V.span(S) True D[3] (3, 0, -5, -6, 0, 0, 1) S.remove(v7) W == V.span(S) True S [(1, 2, 0, -1), (0, -1, 2, 2), (-1, 3, -3, 4)] S == [v1, v3, v4] True Notice that S begins with all seven original vectors, and slowly gets whittled down to just the list [v1, v3, v4]. If you experiment with the above commands, be sure to return to the start and work your way through in order, or the results will not be right. As a bonus, notice that the set of relations of linear dependence provided by Sage, D, is itself a linearly independent set (but within a very different vector space). Is that too weird? U = QQ^7 U.linear_dependence(D) == [] True Now, can you answer the extra credit question from Example COV using Sage? RS Reducing a Span Theorem BS allows us to construct a reduced spanning set for a span. As with the theorem, employing Sage we begin by constructing a matrix with the vectors of the spanning set as columns. Here is a do-over of Example RSC4, illustrating the use of Theorem BS in Sage. V = QQ^4 v1 = vector(QQ, [1,1,2,1]) v2 = vector(QQ, [2,2,4,2]) v3 = vector(QQ, [2,0,-1,1]) v4 = vector(QQ, [7,1,-1,4]) v5 = vector(QQ, [0,2,5,1]) S = [v1, v2, v3, v4, v5] A = column_matrix(S) Chapter V 42 T = [A.column(p) for p in A.pivots()] T [(1, 1, 2, 1), (2, 0, -1, 1)] V.linear_dependence(T) == [] True V.span(S) == V.span(T) True Notice how we compute T with the single line that mirrors the construction of the set T = {vd1 , vd2 , vd3 , . . . vdr } in the statement of Theorem BS. Again, the row- reducing is hidden in the use of the .pivot() matrix method, which necessarily must compute the reduced row-echelon form. The final two compute cells verify both conclusions of the theorem. RES Reworking a Span As another demonstration of using Sage to help us understand spans, linear combi- nations, linear independence and reduced row-echelon form, we will recreate parts of Example RES. Most of this should be familiar, but see the commments following. V = QQ^4 v1 = vector(QQ, [2,1,3,2]) v2 = vector(QQ, [-1,1,0,1]) v3 = vector(QQ, [-8,-1,-9,-4]) v4 = vector(QQ, [3,1,-1,-2]) v5 = vector(QQ, [-10,-1,-1,4]) y = 6*v1 - 7*v2 + v3 +6*v4 + 2*v5 y (9, 2, 1, -3) R = [v1, v2, v3, v4, v5] X = V.span(R) y in X True A = column_matrix(R) P = [A.column(p) for p in A.pivots()] W = V.span(P) W == X True Chapter V 43 y in W True coeff = column_matrix(P) coeff.solve_right(y) (1, -1, 2) coeff.right_kernel() Vector space of degree 3 and dimension 0 over Rational Field Basis matrix: [] V.linear_dependence(P) == [] True The final two results — a trivial null space for coeff and the linear independence of P — both individually imply that the solution to the system of equations (just prior) is unique. Sage produces its own linearly independent spanning set for each span, as we see whenever we inquire about a span. X Vector space of degree 4 and dimension 3 over Rational Field Basis matrix: [ 1 0 0 -8/15] [ 0 1 0 7/15] [ 0 0 1 13/15] Can you extract the three vectors that Sage uses to span X and solve the appro- priate system of equations to see how to write y as a linear combination of these three vectors? Once you have done that, check your answer by hand and think about how using Sage could have been overkill for this question. Section O: Orthogonality EVIC Exact Versus Inexact Computations We are now at a crossroads in our use of Sage. So far our computations have involved rational numbers: fractions of two integers. Sage is able to work with integers of seemingly unlimited size, and then can work with rational numbers exactly. So all of our computations have been exactly correct so far. In practice, many computations, especially those that originate with data, are not so precise. Then we represent real numbers by “floating point numbers.” Since the real numbers are infinite, finite computers must fake it with an extremely large, but still finite, collection of numbers. The price we pay is that some computations will be just slightly imprecise when there is no number available that represents the exact answer. You should now appreciate two problems that occur. If we were to row-reduce a matrix with floating point numbers, there are potentially many computations and Chapter V 44 if a small amount of imprecision arises in each one, these errors can accumulate and lead to wildly incorrect answers. When we row-reduce a matrix, whether or not an entry is zero or not is critically important in the decisions we make about which row operation to perform. If we have an extremely small number (like 10−16 ) how can we be sure if it is zero or not? q 7 Why discuss this now? What is α = 3 ? Hard to say exactly, but it is definitely not a rational number. Norms of vectors will feature prominently in all our discussions about orthogonal vectors, so we now have to recognize the need to work with square roots properly. We have two strategies in Sage. The number system QQbar, also known as the “field of algebraic numbers,” is a truly amazing feature of Sage. It contains the rational numbers, plus every root of every polynomial with coefficients that are rational numbers. For example, notice that α above is one solution to the polynomial equation 3x2 − 7 = 0 and thus is a number in QQbar, so Sage can work with it exactly. These numbers are called “algebraic numbers” and you can recognize them since they print with a question mark near the end to remind you that when printed as a decimal they are approx- imations of numbers that Sage carries internally as exact quantities. For example α can be created with QQbar(sqrt(7/3)) and will print as 1.527525231651947?. Notice that complex numbers begin with the introduction of the imaginary number i, which is a root of the polynomial equation x2 + 1 = 0, so the field of algebraic numbers contains many complex numbers. The downside of QQbar is that com- putations are slow (relatively speaking), so this number system is most useful for examples and demonstrations. The other strategy is to work strictly with approximate numbers, cognizant of the potential for inaccuracies. Sage has two such number systems: RDF and CDF, which are comprised of “double precision” floating point numbers, first limited to just the reals, then expanded to the complexes. Double-precision refers to the use of 64 bits to store the sign, mantissa and exponent in the representation of a real number. This gives 53 bits of precision. Do not confuse these fields with RR and CC, which are similar in appearance but very different in implementation. Sage has implementations of several computations designed exclusively for RDF and CDF, such as the norm. And they are very, very fast. But some computations, like echelon form, can be wildly unreliable with these approximate numbers. We will have more to say about this as we go. In practice, you can use CDF, since RDF is a subset and only different in very limited cases. In summary, QQbar is an extension of QQ which allows exact computations, but can be slow for large examples. RDF and CDF are fast, with special algorithms to control much of the imprecision in some, but not all, computations. So we need to be vigilant and skeptical when we work with these approximate numbers. We will use both strategies, as appropriate. CNIP Conjugates, Norms and Inner Products Conjugates, of complex numbers and of vectors, are straightforward, in QQbar or in CDF. alpha = QQbar(2 + 3*I) alpha.conjugate() 2 - 3*I beta = CDF(2+3*I) beta.conjugate() 2.0 - 3.0*I Chapter V 45 v = vector(QQbar, [5-3*I, 2+6*I]) v.conjugate() (5 + 3*I, 2 - 6*I) w = vector(CDF, [5-3*I, 2+6*I]) w.conjugate() (5.0 + 3.0*I, 2.0 - 6.0*I) The term “inner product” means slightly different things to different people. For some, it is the “dot product” that you may have seen in a calculus or physics course. Our inner product could be called the “Hermitian inner product” to emphasize the use of vectors over the complex numbers and conjugating some of the entries. So Sage has a .dot_product(), .inner_product(), and .hermitian_inner_product() — we want to use the last one. From now on, when we mention an inner product in the context of using Sage, we will mean .hermitian_inner_product(). We will redo the first part of Example CSIP. Notice that the syntax is a bit asymmetric. u = vector(QQbar, [2+3*I, 5+2*I, -3+I]) v = vector(QQbar, [1+2*I, -4+5*I, 5*I]) u.hermitian_inner_product(v) 3 + 19*I Norms are as easy as conjugates. Easier maybe. It might be useful to realize that Sage uses entirely distinct code to compute an exact norm over QQbar versus an approximate norm over CDF, though that is totally transparent as you issue commands. Here is Example CNSV reprised. entries = [3+2*I, 1-6*I, 2+4*I, 2+I] u = vector(QQbar, entries) u.norm() 8.66025403784439? u = vector(CDF, entries) u.norm() 8.66025403784 numerical_approx(5*sqrt(3), digits = 30) 8.66025403784438646763723170753 We have three different numerical approximations, the latter 30-digit number being an approximation to the answer in the text. But there is no inconsistency between them. The first, an algebraic number, is represented internally as 5 ∗ a where √ a is a root of the polynomial equation x2 − 3 = 0, in other words it is 5 3. The CDF value prints with a few digits less than what is carried internally. Notice that our different definitions of the inner product make no difference in the computation of a norm. Chapter V 46 One warning now that we are working with complex numbers. It is easy to “clobber” the symbol I used for the imaginary number i. In other words, Sage will allow you to assign it to something else, rendering it useless. An identity matrix is a likely reassignment. If you run the next compute cell, be sure to evaluate the compute cell afterward to restore I to its usual role. alpha = QQbar(5 - 6*I) I = identity_matrix(2) beta = QQbar(2+5*I) Traceback (most recent call last): ... TypeError: Illegal initializer for algebraic number restore() I^2 -1 We will finish with a verification of Theorem IPN. To test equality it is best if we work with entries from QQbar. v = vector(QQbar, [2-3*I, 9+5*I, 6+2*I, 4-7*I]) v.hermitian_inner_product(v) == v.norm()^2 True OGS Orthogonality and Gram-Schmidt It is easy enough to check a pair of vectors for orthogonality (is the inner product zero?). To check that a set is orthogonal, we just need to do this repeatedly. This is a redo of Example AOS. x1 = vector(QQbar, [ 1+I, 1, 1-I, I]) x2 = vector(QQbar, [ 1+5*I, 6+5*I, -7-I, 1-6*I]) x3 = vector(QQbar, [-7+34*I, -8-23*I, -10+22*I, 30+13*I]) x4 = vector(QQbar, [ -2-4*I, 6+I, 4+3*I, 6-I]) S = [x1, x2, x3, x4] ips = [S[i].hermitian_inner_product(S[j]) for i in range(3) for j in range(i+1,3)] all([ip == 0 for ip in ips]) True Notice how the list comprehension computes each pair just once, and never checks the inner product of a vector with itself. If we wanted to check that a set is orthonormal, the “normal” part is less involved. We will check the set above, even though we can clearly see that the four vectors are not even close to being unit vectors. Be sure to run the above definitions of S before running the next compute cell. ips = [S[i].hermitian_inner_product(S[i]) for i in range(3)] all([ip == 1 for ip in ips]) Chapter V 47 False Applying the Gram-Schmidt procedure to a set of vectors is the type of compu- tation that a program like Sage is perfect for. Gram-Schmidt is implemented as a method for matrices, where we interpret the rows of the matrix as the vectors in the original set. The result is two matrices, where the first has rows that are the orthogonal vectors. The second matrix has rows that provide linear combinations of the orthogonal vectors that equal the original vectors. The original vectors do not need to form a linearly independent set, and when the set is linearly dependent, then zero vectors produced are not part of the returned set. Over CDF the set is automatically orthonormal, and since a different algorithm is used (to help control the imprecisions), the results will look different than what would result from Theorem GSP. We will illustrate with the vectors from Example GSTV. v1 = vector(CDF, [ 1, 1+I, 1]) v2 = vector(CDF, [-I, 1, 1+I]) v3 = vector(CDF, [ 0, I, I]) A = matrix([v1,v2,v3]) G, M = A.gram_schmidt() G.round(5) [ -0.5 -0.5 - 0.5*I -0.5] [ 0.30151 + 0.45227*I -0.15076 + 0.15076*I -0.30151 - 0.75378*I] [ 0.6396 + 0.2132*I -0.2132 - 0.6396*I 0.2132 + 0.2132*I] We formed the matrix A with the three vectors as rows, and of the two outputs we are interested in the first one, whose rows form the orthonormal set. We round the numbers to 5 digits, just to make the result fit nicely on your screen. Let us do it again, now exactly over QQbar. We will output the entries of the matrix as list, working across rows first, so it fits nicely. v1 = vector(QQbar, [ 1, 1+I, 1]) v2 = vector(QQbar, [-I, 1, 1+I]) v3 = vector(QQbar, [ 0, I, I]) A = matrix([v1,v2,v3]) G, M = A.gram_schmidt(orthonormal=True) Sequence(G.list(), cr=True) [ 0.50000000000000000?, 0.50000000000000000? + 0.50000000000000000?*I, 0.50000000000000000?, -0.3015113445777636? - 0.4522670168666454?*I, 0.1507556722888818? - 0.1507556722888818?*I, 0.3015113445777636? + 0.7537783614444091?*I, -0.6396021490668313? - 0.2132007163556105?*I, 0.2132007163556105? + 0.6396021490668313?*I, -0.2132007163556105? - 0.2132007163556105?*I ] Notice that we asked for orthonormal output, so the rows of G are the vectors {w1 , w2 , w3 } in Example ONTV. Exactly. We can restrict ourselves to QQ and forego the “normality” to obtain just the orthogonal set {u1 , u2 , u3 } of Example GSTV. Chapter V 48 v1 = vector(QQbar, [ 1, 1+I, 1]) v2 = vector(QQbar, [-I, 1, 1+I]) v3 = vector(QQbar, [ 0, I, I]) A = matrix([v1, v2, v3]) G, M = A.gram_schmidt(orthonormal=False) Sequence(G.list(), cr=True) [ 1, I + 1, 1, -0.50000000000000000? - 0.75000000000000000?*I, 0.25000000000000000? - 0.25000000000000000?*I, 0.50000000000000000? + 1.2500000000000000?*I, -0.2727272727272728? - 0.0909090909090909?*I, 0.0909090909090909? + 0.2727272727272728?*I, -0.0909090909090909? - 0.0909090909090909?*I ] Notice that it is an error to ask for an orthonormal set over QQ since you cannot expect to take square roots of rationals and stick with rationals. v1 = vector(QQ, [1, 1]) v2 = vector(QQ, [2, 3]) A = matrix([v1,v2]) G, M = A.gram_schmidt(orthonormal=True) Traceback (most recent call last): ... TypeError: QR decomposition unable to compute square roots in Rational Field Chapter M Matrices Section MO: Matrix Operations MS Matrix Spaces Sage defines our set Mmn as a “matrix space” with the command MatrixSpace(R, m, n) where R is a number system and m and n are the number of rows and number of columns, respectively. This object does not have much functionality defined in Sage. Its main purposes are to provide a parent for matrices, and to provide another way to create matrices. The two matrix operations just defined (addition and scalar multiplication) are implemented as you would expect. In the example below, we create two matrices in M23 from just a list of 6 entries, by coercing the list into a matrix by using the relevant matrix space as if it were a function. Then we can perform the basic operations of matrix addition (Definition MA) and matrix scalar multiplication (Definition MSM). MS = MatrixSpace(QQ, 2, 3) MS Full MatrixSpace of 2 by 3 dense matrices over Rational Field A = MS([1, 2, 1, 4, 5, 4]) B = MS([1/1, 1/2, 1/3, 1/4, 1/5, 1/6]) A + B [ 2 5/2 4/3] [17/4 26/5 25/6] 60*B [60 30 20] [15 12 10] 5*A - 120*B [-115 -50 -35] [ -10 1 0] Coercion can make some interesting conveniences possible. Notice how the scalar 37 in the following expression is coerced to 37 times an identity matrix of the proper size. 49 Chapter M 50 A = matrix(QQ, [[ 0, 2, 4], [ 6, 0, 8], [10, 12, 0]]) A + 37 [37 2 4] [ 6 37 8] [10 12 37] This coercion only applies to sums with square matrices. You might try this again, but with a rectangular matrix, just to see what the error message says. MO Matrix Operations Every operation in this section is implemented in Sage. The only real subtlety is determining if certain matrices are symmetric, which we will discuss below. In linear algebra, the term “adjoint” has two unrelated meanings, so you need to be careful when you see this term. In particular, in Sage it is used to mean some- thing different. So our version of the adjoint is implemented as the matrix method .conjugate_transpose(). Here are some straightforward examples. A = matrix(QQ, [[-1, 2, 4], [ 0, 3, 1]]) A [-1 2 4] [ 0 3 1] A.transpose() [-1 0] [ 2 3] [ 4 1] A.is_symmetric() False B = matrix(QQ, [[ 1, 2, -1], [ 2, 3, 4], [-1, 4, -6]]) B.is_symmetric() True C = matrix(QQbar, [[ 2-I, 3+4*I], [5+2*I, 6]]) C.conjugate() Chapter M 51 [2 + 1*I 3 - 4*I] [5 - 2*I 6] C.conjugate_transpose() [2 + 1*I 5 - 2*I] [3 - 4*I 6] With these constructions, we can test, or demonstrate, some of the theorems above. Of course, this does not make the theorems true, but is satisfying nonethe- less. This can be an effective technique when you are learning new Sage commands or new linear algebra — if your computations are not consistent with theorems, then your understanding of the linear algebra may be flawed, or your understand- ing of Sage may be flawed, or Sage may have a bug! Note in the following how we use comparison (==) between matrices as an implementation of matrix equality (Definition ME). A = matrix(QQ, [[ 1, -1, 3], [-3, 2, 0]]) B = matrix(QQ, [[5, -2, 7], [1, 3, -2]]) C = matrix(QQbar, [[2+3*I, 1 - 6*I], [3, 5+2*I]]) A.transpose().transpose() == A True (A+B).transpose() == A.transpose() + B.transpose() True (2*C).conjugate() == 2*C.conjugate() True a = QQbar(3 + 4*I) acon = a.conjugate() (a*C).conjugate_transpose() == acon*C.conjugate_transpose() True The opposite is true — you can use theorems to convert, or express, Sage code into alternative, but mathematically equivalent forms. Here is the subtlety. With approximate numbers, such as in RDF and CDF, it can be tricky to decide if two numbers are equal, or if a very small number is zero or not. In these situations Sage allows us to specify a “tolerance” — the largest number that can be effectively considered zero. Consider the following: A = matrix(CDF, [[1.0, 0.0], [0.0, 1.0]]) A Chapter M 52 [1.0 0.0] [0.0 1.0] A.is_symmetric() True A[0,1] = 0.000000000002 A [ 1.0 2e-12] [ 0.0 1.0] A.is_symmetric() False A[0,1] = 0.000000000001 A [ 1.0 1e-12] [ 0.0 1.0] A.is_symmetric() True Clearly the last result is not correct. This is because 0.000000000001 = 1.0 × 10−12 is “small enough” to be confused as equal to the zero in the other corner of the matrix. However, Sage will let us set our own idea of when two numbers are equal, by setting a tolerance on the difference between two numbers that will allow them to be considered equal. The default tolerance is set at 1.0 × 10−12 . Here we use Sage’s syntax for scientific notation to specify the tolerance. A = matrix(CDF, [[1.0, 0.0], [0.0, 1.0]]) A.is_symmetric() True A[0,1] = 0.000000000001 A.is_symmetric() True A.is_symmetric(tol=1.0e-13) Chapter M 53 False This is not a course in numerical linear algebra, even if that is a fascinating field of study. To concentrate on the main ideas of introductory linear algebra, whenever possible we will concentrate on number systems like the rational numbers or algebraic numbers where we can rely on exact results. If you are ever unsure if a number system is exact or not, just ask. QQ.is_exact() True RDF.is_exact() False Section MM: Matrix Multiplication MVP Matrix-Vector Product A matrix-vector product is very natural in Sage, and we can check the result against a linear combination of the columns. A = matrix(QQ, [[1, -3, 4, 5], [2, 3, -2, 0], [5, 6, 8, -2]]) v = vector(QQ, [2, -2, 1, 3]) A*v (27, -4, 0) sum([v[i]*A.column(i) for i in range(len(v))]) (27, -4, 0) Notice that when a matrix is square, a vector of the correct size can be used in Sage in a product with a matrix by placing the vector on either side of the matrix. However, the two results are not the same, and we will not have ocassion to place the vector on the left any time soon. So, despite the possibility, be certain to keep your vectors on the right side of a matrix in a product. B = matrix(QQ, [[ 1, -3, 4, 5], [ 2, 3, -2, 0], [ 5, 6, 8, -2], [-4, 1, 1, 2]]) w = vector(QQ, [1, 2, -3, 2]) B*w (-7, 14, -11, -1) w*B Chapter M 54 (-18, -13, -22, 15) B*w == w*B False Since a matrix-vector product forms a linear combination of columns of a matrix, it is now very easy to check if a vector is a solution to a system of equations. This is basically the substance of Theorem SLEMM. Here we construct a system of equations and construct two solutions and one non-solution by applying Theorem PSPHS. Then we use a matrix-vector product to verify that the vectors are, or are not, solutions. coeff = matrix(QQ, [[-1, 3, -1, -1, 0, 2], [ 2, -6, 1, -2, -5, -8], [ 1, -3, 2, 5, 4, 1], [ 2, -6, 2, 2, 1, -3]]) const = vector(QQ, [13, -25, -17, -23]) solution1 = coeff.solve_right(const) coeff*solution1 (13, -25, -17, -23) nsp = coeff.right_kernel(basis=’pivot’) nsp Vector space of degree 6 and dimension 3 over Rational Field User basis matrix: [ 3 1 0 0 0 0] [ 3 0 -4 1 0 0] [ 1 0 1 0 -1 1] nspb = nsp.basis() solution2 = solution1 + 5*nspb[0]+(-4)*nspb[1]+2*nspb[2] coeff*solution2 (13, -25, -17, -23) nonnullspace = vector(QQ, [5, 0, 0, 0, 0, 0]) nonnullspace in nsp False nonsolution = solution1 + nonnullspace coeff*nonsolution (8, -15, -12, -13) We can now explain the difference between “left” and “right” variants of various Sage commands for linear algebra. Generally, the direction refers to where the vector is placed in a matrix-vector product. We place a vector on the right and understand Chapter M 55 this to mean a linear combination of the columns of the matrix. Placing a vector to the left of a matrix can be understood, in a manner totally consistent with our upcoming definition of matrix multiplication, as a linear combination of the rows of the matrix. So the difference between A.solve_right(v) and A.solve_left(v) is that the former asks for a vector x such that A*x == v, while the latter asks for a vec- tor x such that x*A == v. Given Sage’s preference for rows, a direction-neutral version of a command, if it exists, will be the “left” version. For example, there is a .right_kernel() matrix method, while the .left_kernel() and .kernel() methods are identical — the names are synonyms for the exact same routine. So when you see a Sage command that comes in “left” and “right” variants figure out just what part of the defined object involves a matrix-vector product and form an interpretation from that. MM Matrix Multiplication Matrix multiplication is very natural in Sage, and is just as easy as multiplying two numbers. We illustrate Theorem EMP by using it to compute the entry in the first row and third column. A = matrix(QQ, [[3, -1, 2, 5], [9, 1, 2, -4]]) B = matrix(QQ, [[1, 6, 1], [0, -1, 2], [5, 2, 3], [1, 1, 1]]) A*B [18 28 12] [15 53 13] sum([A[0,k]*B[k,2] for k in range(A.ncols())]) 12 Note in the final statement, we could replace A.ncols() by B.nrows() since these two quantities must be identical. You can experiment with the last statement by editing it to compute any of the five other entries of the matrix product. Square matrices can be multiplied in either order, but the result will almost always be different. Execute repeatedly the following products of two random 4 × 4 matrices, with a check on the equality of the two products in either order. It is possible, but highly unlikely, that the two products will be equal. So if this compute cell ever produces True it will be a minor miracle. A = random_matrix(QQ,4,4) B = random_matrix(QQ,4,4) A*B == B*A # random, sort of False PMM Properties of Matrix Multiplication As before, we can use Sage to demonstrate theorems. With randomly-generated matrices, these verifications might be even more believable. Some of the above results should feel fairly routine, but some are perhaps contrary to intuition. For example, Theorem MMT might at first glance seem surprising due to the require- ment that the order of the product is reversed. Here is how we would investigate this theorem in Sage. The following compute cell should always return True. Repeated experimental evidence does not make a proof, but certainly gives us confidence. Chapter M 56 A = random_matrix(QQ, 3, 7) B = random_matrix(QQ, 7, 5) (A*B).transpose() == B.transpose()*A.transpose() True By now, you can probably guess the matrix method for checking if a matrix is Hermitian. A = matrix(QQbar, [[ 45, -5-12*I, -1-15*I, -56-8*I], [-5+12*I, 42, 32*I, -14-8*I], [-1+15*I, -32*I, 57, 12+I], [-56+8*I, -14+8*I, 12-I, 93]]) A.is_hermitian() True We can illustrate the most fundamental property of a Hermitian matrix. The vectors x and y below are random, but according to Theorem HMIP the final command should produce True for any possible values of these two vectors. (You would be right to think that using random vectors over QQbar would be a better idea, but at this writing, these vectors are not as “random” as one would like, and are insufficient to perform an accurate test here.) x = random_vector(QQ, 4) + QQbar(I)*random_vector(QQ, 4) y = random_vector(QQ, 4) + QQbar(I)*random_vector(QQ, 4) (A*x).hermitian_inner_product(y) == x.hermitian_inner_product(A*y) True Section MISLE: Matrix Inverses and Systems of Linear Equations MISLE Matrix Inverse, Systems of Equations We can use the computational method described in this section in hopes of finding a matrix inverse, as Theorem CINM gets us halfway there. We will continue with the matrix from Example MI. First we check that the matrix is nonsingular so we can apply the theorem, then we get “half” an inverse, and verify that it also behaves as a “full” inverse by meeting the full definition of a matrix inverse (Definition MI). A = matrix(QQ, [[ 1, 2, 1, 2, 1], [-2, -3, 0, -5, -1], [ 1, 1, 0, 2, 1], [-2, -3, -1, -3, -2], [-1, -3, -1, -3, 1]]) A.is_singular() False I5 = identity_matrix(5) M = A.augment(I5); M Chapter M 57 [ 1 2 1 2 1 1 0 0 0 0] [-2 -3 0 -5 -1 0 1 0 0 0] [ 1 1 0 2 1 0 0 1 0 0] [-2 -3 -1 -3 -2 0 0 0 1 0] [-1 -3 -1 -3 1 0 0 0 0 1] N = M.rref(); N [ 1 0 0 0 0 -3 3 6 -1 -2] [ 0 1 0 0 0 0 -2 -5 -1 1] [ 0 0 1 0 0 1 2 4 1 -1] [ 0 0 0 1 0 1 0 1 1 0] [ 0 0 0 0 1 1 -1 -2 0 1] J = N.matrix_from_columns(range(5,10)); J [-3 3 6 -1 -2] [ 0 -2 -5 -1 1] [ 1 2 4 1 -1] [ 1 0 1 1 0] [ 1 -1 -2 0 1] A*J == I5 True J*A == I5 True Note that the matrix J is constructed by taking the last 5 columns of N (num- bered 5 through 9) and using them in the matrix_from_columns() matrix method. What happens if you apply the procedure above to a singular matrix? That would be an instructive experiment to conduct. With an inverse of a coefficient matrix in hand, we can easily solve systems of equations, in the style of Example SABMI. We will recycle the matrices A and its inverse, J, from above, so be sure to run those compute cells first if you are playing along. We consider a system with A as a coefficient matrix and solve a linear system twice, once the old way and once the new way. Recall that with a nonsingular coefficient matrix, the solution will be unique for any choice of const, so you can experiment by changing the vector of constants and re-executing the code. const = vector(QQ, [3, -4, 2, 1, 1]) A.solve_right(const) (-12, -2, 3, 6, 4) J*const Chapter M 58 (-12, -2, 3, 6, 4) A.solve_right(const) == J*const True Section MINM: Matrix Inverses and Nonsingular Matrices MI Matrix Inverse Now we know that invertibility is equivalent to nonsingularity, and that the pro- cedure outlined in Theorem CINM will always yield an inverse for a nonsingular matrix. But rather than using that procedure, Sage implements a .inverse() method. In the following, we compute the inverse of a 3 × 3 matrix, and then purposely convert it to a singular matrix by replacing the last column by a linear combination of the first two. A = matrix(QQ,[[ 3, 7, -6], [ 2, 5, -5], [-3, -8, 8]]) A.is_singular() False Ainv = A.inverse(); Ainv [ 0 8 5] [ 1 -6 -3] [ 1 -3 -1] A*Ainv [1 0 0] [0 1 0] [0 0 1] col_0 = A.column(0) col_1 = A.column(1) C = column_matrix([col_0, col_1, 2*col_0 - 4*col_1]) C.is_singular() True C.inverse() Traceback (most recent call last): ... ZeroDivisionError: input matrix must be nonsingular Chapter M 59 Notice how the failure to invert C is explained by the matrix being singular. Systems with nonsingular coefficient matrices can be solved easily with the ma- trix inverse. We will recycle A as a coefficient matrix, so be sure to execute the code above. const = vector(QQ, [2, -1, 4]) A.solve_right(const) (12, -4, 1) A.inverse()*const (12, -4, 1) A.solve_right(const) == A.inverse()*const True If you find it more convenient, you can use the same notation as the text for a matrix inverse. A^-1 [ 0 8 5] [ 1 -6 -3] [ 1 -3 -1] NME3 Nonsingular Matrix Equivalences, Round 3 For square matrices, Sage has the methods .is_singular() and .is_invertible(). By Theorem NI we know these two functions to be logical opposites. One way to express this is that these two methods will always return different values. Here we demonstrate with a nonsingular matrix and a singular matrix. The comparison =! is “not equal.” nonsing = matrix(QQ, [[ 0, -1, 1, -3], [ 1, 1, 0, -3], [ 0, 4, -3, 8], [-2, -4, 1, 5]]) nonsing.is_singular() != nonsing.is_invertible() True sing = matrix(QQ, [[ 1, -2, -6, 8], [-1, 3, 7, -8], [ 0, -4, -3, -2], [ 0, 3, 1, 4]]) sing.is_singular() != sing.is_invertible() True We could test other properties of the matrix inverse, such as Theorem SS. Chapter M 60 A = matrix(QQ, [[ 3, -5, -2, 8], [-1, 2, 0, -1], [-2, 4, 1, -4], [ 4, -5, 0, 8]]) B = matrix(QQ, [[ 1, 2, 4, -1], [-2, -3, -8, -2], [-2, -4, -7, 5], [ 2, 5, 7, -8]]) A.is_invertible() and B.is_invertible() True (A*B).inverse() == B.inverse()*A.inverse() True UM Unitary Matrices No surprise about how we check if a matrix is unitary. Here is Example UM3, A = matrix(QQbar, [ [(1+I)/sqrt(5), (3+2*I)/sqrt(55), (2+2*I)/sqrt(22)], [(1-I)/sqrt(5), (2+2*I)/sqrt(55), (-3+I)/sqrt(22)], [ I/sqrt(5), (3-5*I)/sqrt(55), (-2)/sqrt(22)] ]) A.is_unitary() True A.conjugate_transpose() == A.inverse() True We can verify Theorem UMPIP, where the vectors u and v are created randomly. Try evaluating this compute cell with your own choices. u = random_vector(QQ, 3) + QQbar(I)*random_vector(QQ, 3) v = random_vector(QQ, 3) + QQbar(I)*random_vector(QQ, 3) (A*u).hermitian_inner_product(A*v) == u.hermitian_inner_product(v) True If you want to experiment with permutation matrices, Sage has these too. We can create a permutation matrix from a list that indicates for each column the row with a one in it. Notice that the product here of two permutation matrices is again a permutation matrix. sigma = Permutation([2,3,1]) S = sigma.to_matrix(); S [0 0 1] [1 0 0] [0 1 0] Chapter M 61 tau = Permutation([1,3,2]) T = tau.to_matrix(); T [1 0 0] [0 0 1] [0 1 0] S*T [0 1 0] [1 0 0] [0 0 1] rho = Permutation([2, 1, 3]) S*T == rho.to_matrix() True Section CRS: Column and Row Spaces CSCS Column Space and Consistent Systems We could compute the column space of a matrix with a span of the set of columns of the matrix, much as we did back in Sage CSS when we were checking consis- tency of linear systems using spans of the set of columns of a coefficent matrix. However, Sage provides a convenient matrix method to construct this same span: .column_space(). Here is a check. D = matrix(QQ, [[ 2, -1, -4], [ 5, 2, -1], [-3, 1, 5]]) cs = D.column_space(); cs Vector space of degree 3 and dimension 2 over Rational Field Basis matrix: [ 1 0 -11/9] [ 0 1 -1/9] cs_span = (QQ^3).span(D.columns()) cs == cs_span True In Sage CSS we discussed solutions to systems of equations and the membership of the vector of constants in the span of the columns of the coefficient matrix. This discussion turned on Theorem SLSLC. We can now be slightly more efficient with Theorem CSCS, while still using the same ideas. We recycle the computations from Example CSMCS that use Archetype D and Archetype E. Chapter M 62 coeff = matrix(QQ, [[ 2, 1, 7, -7], [-3, 4, -5, -6], [ 1, 1, 4, -5]]) constD = vector(QQ, [8, -12, 4]) constE = vector(QQ, [2, 3, 2]) cs = coeff.column_space() coeff.solve_right(constD) (4, 0, 0, 0) constD in cs True coeff.solve_right(constE) Traceback (most recent call last): ... ValueError: matrix equation has no solutions constE in cs False CSOC Column Space, Original Columns It might be worthwhile for Sage to create a column space using actual columns of the matrix as a spanning set. But we can do it ourselves fairly easily. A discussion follows the example. A = matrix(QQ, [[-1, -1, 0, 1, 0, 1, -2, -3, 0], [-1, 0, 0, 5, -2, 6, -6, 3, 5], [ 0, 0, 1, -6, -1, 0, -5, 0, 5], [ 2, 2, 1, -8, -2, 0, -4, 8, 8], [ 2, 2, 0, -2, -1, 0, 1, 8, 3], [ 2, 1, 0, -6, -1, -1, -1, 6, 4]]) A.rref() [ 1 0 0 -5 0 -2 0 1 1] [ 0 1 0 4 0 1 2 2 -1] [ 0 0 1 -6 0 -2 -2 -2 2] [ 0 0 0 0 1 -2 3 -2 -3] [ 0 0 0 0 0 0 0 0 0] [ 0 0 0 0 0 0 0 0 0] B = A.matrix_from_columns(A.pivots()) A.column_space() == B.column_space() True Chapter M 63 B.rref() [1 0 0 0] [0 1 0 0] [0 0 1 0] [0 0 0 1] [0 0 0 0] [0 0 0 0] cs = (QQ^6).span_of_basis(B.columns()) cs Vector space of degree 6 and dimension 4 over Rational Field User basis matrix: [-1 -1 0 2 2 2] [-1 0 0 2 2 1] [ 0 0 1 1 0 0] [ 0 -2 -1 -2 -1 -1] cs.basis() [ (-1, -1, 0, 2, 2, 2), (-1, 0, 0, 2, 2, 1), (0, 0, 1, 1, 0, 0), (0, -2, -1, -2, -1, -1) ] cs.echelonized_basis() [ (1, 0, 0, 0, 0, 5), (0, 1, 0, 0, 0, -1), (0, 0, 1, 0, -1, -3), (0, 0, 0, 1, 1, 3) ] cs == A.column_space() True cs2 = (QQ^6).span_of_basis([A.column(i) for i in A.pivots()]) cs2 == A.column_space() True We see that A has four pivot columns, numbered 0,1,2,4. The matrix B is just a convenience to hold the pivot columns of A. However, the column spaces of A and B should be equal, as Sage verifies. Also B will row-reduce to the same 0-1 pivot columns of the reduced row-echelon form of the full matrix A. So it is no Chapter M 64 accident that the reduced row-echelon form of B is a full identity matrix, followed by sufficiently many zero rows to give the matrix the correct size. The vector space method .span_of_basis() is new to us. It creates a span of a set of vectors, as before, but we now are responsible for supplying a linearly independent set of vectors. Which we have done. We know this because Theorem BCS guarantees the set we provided is linearly independent (and spans the column space), while Sage would have given us an error if we had provided a linearly dependent set. In return, Sage will carry this linearly independent spanning set along with the vector space, something Sage calls a “user basis.” Notice how cs has two linearly independent spanning sets now. Our set of “origi- nal columns” is obtained via the standard vector space method .basis() and we can obtain a linearly independent spanning set that looks more familiar with the vector space method .echelonized_basis(). For a vector space created with a simple .span() construction these two commands would yield identical results — it is only when we supply a linearly independent spanning set with the .span_of_basis() method that a “user basis” becomes relevant. Finally, we check that cs is indeed the column space of A (we knew it would be) and then we provide a one-line, totally general construction of the column space using original columns. This is an opportunity to make an interesting observation, which could be used to substantiate several theorems. When we take the original columns that we rec- ognize as pivot columns, and use them alone to form a matrix, this new matrix will always row-reduce to an identity matrix followed by zero rows. This is basically a consequence of reduced row-echelon form. Evaluate the compute cell below repeat- edly. The number of columns could in theory change, though this is unlikely since the columns of a random matrix are unlikely to be linearly dependent. In any event, the form of the result will always be an identity matrix followed by some zero rows. F = random_matrix(QQ, 5, 3) F.matrix_from_columns(F.pivots()).rref() # random [1 0 0] [0 1 0] [0 0 1] [0 0 0] [0 0 0] With more columns than rows, we know by Theorem MVSLD that we will have a reduced number of pivot columns. Here, we will almost always see an identity matrix as the result, though we could get a smaller identity matrix followed by zero rows. F = random_matrix(QQ, 3, 5) F.matrix_from_columns(F.pivots()).rref() # random [1 0 0] [0 1 0] [0 0 1] NME4 Nonsingular Matrix Equivalences, Round 4 Archetype A and Archetype B have square coefficient matrices that illustrate the di- chotomy of singular and nonsingular matrices. Here we illustrate the latest addition to our series of equivalences, Theorem CSNM. A = matrix(QQ, [[1, -1, 2], [2, 1, 1], [1, 1, 0]]) Chapter M 65 B = matrix(QQ, [[-7, -6, -12], [ 5, 5, 7], [ 1, 0, 4]]) A.is_singular() True A.column_space() == QQ^3 False B.is_singular() False B.column_space() == QQ^3 True We can even write Theorem CSNM as a one-line Sage statement that will return True for any square matrix. Run the following repeatedly and it should always return True. We have kept the size of the matrix relatively small to be sure that some of the random matrices produced are singular. A = random_matrix(QQ, 4, 4) A.is_singular() == (not A.column_space() == QQ^4) True RSM Row Space of a Matrix Not to be outdone, and not suprisingly, Sage can compute a row space with the matrix method .row_space(). Indeed, given Sage’s penchant for treating vectors as rows, much of Sage’s infrastructure for vector spaces ultimately relies on Theorem REMRS. More on that in Sage SUTH0. For now, we reprise Example IAS as an illustration. v1 = vector(QQ, [1, 2, 1, 6, 6]) v2 = vector(QQ, [3, -1, 2, -1, 6]) v3 = vector(QQ, [1, -1, 0, -1, -2]) v4 = vector(QQ, [-3, 2, -3, 6, -10]) X = (QQ^5).span([v1, v2, v3, v4]) A = matrix([v1, v2, v3, v4]) rsA = A.row_space() X == rsA True B = A.rref() rsB = B.row_space() X == rsB Chapter M 66 True X Vector space of degree 5 and dimension 3 over Rational Field Basis matrix: [ 1 0 0 2 -1] [ 0 1 0 3 1] [ 0 0 1 -2 5] X.basis() [ (1, 0, 0, 2, -1), (0, 1, 0, 3, 1), (0, 0, 1, -2, 5) ] B [ 1 0 0 2 -1] [ 0 1 0 3 1] [ 0 0 1 -2 5] [ 0 0 0 0 0] We begin with the same four vectors in Example IAS and create their span, the vector space X. The matrix A has these four vectors as rows and B is the reduced row-echelon form of A. Then the row spaces of A and B are equal to the vector space X (and each other). The way Sage describes this vector space is with a matrix whose rows are the nonzero rows of the reduced row-echelon form of the matrix A. This is Theorem BRS in action where we go with Sage’s penchant for rows and ignore the text’s penchant for columns. We can illustrate a few other results about row spaces with Sage. Discussion follows. A = matrix(QQ,[[1, 1, 0, -1, 1, 1, 0], [4, 5, 1, -6, 1, 6, 1], [5, 5, 1, -5, 4, 5, 2], [3, -1, 0, 5, 11, -5, 4]]) A.row_space() == A.transpose().column_space() True B = matrix(QQ,[[ 7, 9, 2, -11, 1, 11, 2], [-4, -3, 1, 2, -7, -2, 1], [16, 8, 2, 0, 30, 0, 12], [ 2, 10, 2, -18, -16, 18, -4]]) B.column_space() == B.transpose().row_space() Chapter M 67 True A.rref() == B.rref() True A.row_space() == B.row_space() True We use the matrix A to illustrate Definition RSM, and the matrix B to illustrate Theorem CSRST. A and B were designed to have the same reduced row-echelon form, and hence be row-equivalent, so this is not a consequence of any theorem or previous computation. However, then Theorem REMRS guarantees that the row spaces of A and B are equal. Section FS: Four Subsets LNS Left Null Spaces Similar to (right) null spaces, a left null space can be computed in Sage with the matrix method .left_kernel(). For a matrix A, elements of the (right) null space are vectors x such that Ax = 0. If we interpret a vector placed to the left of a matrix as a 1-row matrix, then the product xA can be interpreted, according to our definition of matrix multiplication (Definition MM, as the entries of the vector x providing the scalars for a linear combination of the rows of the matrix A. So you can view each vector in the left null space naturally as the entries of a matrix with a single row, Y , with the property that Y A = O. Given Sage’s preference for row vectors, the simpler name .kernel() is a syn- onym for .left_kernel(). Given your text’s preference for column vectors, we will continue to rely on the .right_kernel(). Left kernels in Sage have the same options as right kernels and produce vector spaces as output in the same manner. Once created, a vector space all by itself (with no history) is neither left or right. Here is a quick repeat of Example LNS. A = matrix(QQ, [[ 1, -3, 1], [-2, 1, 1], [ 1, 5, 1], [ 9, -4, 0]]) A.left_kernel(basis=’pivot’) Vector space of degree 4 and dimension 1 over Rational Field User basis matrix: [-2 3 -1 1] RRSM Row-Reducing a Symbolic Matrix Sage can very nearly reproduce the reduced row-echelon form we obtained from the augmented matrix with variables present in the final column. The first line below is a bit of advanced Sage. It creates a number system R mixing the rational numbers with the variables b1 through b6. It is not important to know the details of this construction right now. B is the reduced row-echelon form of A, as computed by Sage, where we have displayed the last column separately so it will all fit. You will notice that B is different than what we used in Example CSANS, where all the differences are in the final column. Chapter M 68 However, we can perform row operations on the final two rows of B to bring Sage’s result in line with what we used above for the final two entries of the last column, which are the most critical. Notice that since the final two rows are almost all zeros, any sequence of row operations on just these two rows will preserve the zeros (and we need only display the final column to keep track of our progress). R.<b1,b2,b3,b4,b5,b6> = QQ[] A = matrix(R, [[ 10, 0, 3, 8, 7, b1], [-16, -1, -4, -10, -13, b2], [ -6, 1, -3, -6, -6, b3], [ 0, 2, -2, -3, -2, b4], [ 3, 0, 1, 2, 3, b5], [ -1, -1, 1, 1, 0, b6]]) A [ 10 0 3 8 7 b1] [-16 -1 -4 -10 -13 b2] [ -6 1 -3 -6 -6 b3] [ 0 2 -2 -3 -2 b4] [ 3 0 1 2 3 b5] [ -1 -1 1 1 0 b6] B = copy(A.echelon_form()) B[0:6,0:5] [ 1 0 0 0 2] [ 0 1 0 0 -3] [ 0 0 1 0 1] [ 0 0 0 1 -2] [ 0 0 0 0 0] [ 0 0 0 0 0] B[0:6,5] [ -1/4*b1 - 5/4*b2 + 11/4*b3 - 2*b4] [ 3*b2 - 8*b3 + 6*b4] [ -3/2*b1 + 3/2*b2 - 13/2*b3 + 4*b4] [ b1 + b2 - b3 + b4] [ 1/4*b1 + 1/4*b2 + 1/4*b3 + b5] [1/4*b1 - 3/4*b2 + 9/4*b3 - b4 + b6] B.rescale_row(4, 4); B[0:6, 5] [ -1/4*b1 - 5/4*b2 + 11/4*b3 - 2*b4] [ 3*b2 - 8*b3 + 6*b4] [ -3/2*b1 + 3/2*b2 - 13/2*b3 + 4*b4] [ b1 + b2 - b3 + b4] [ b1 + b2 + b3 + 4*b5] [1/4*b1 - 3/4*b2 + 9/4*b3 - b4 + b6] Chapter M 69 B.add_multiple_of_row(5, 4, -1/4); B[0:6, 5] [-1/4*b1 - 5/4*b2 + 11/4*b3 - 2*b4] [ 3*b2 - 8*b3 + 6*b4] [-3/2*b1 + 3/2*b2 - 13/2*b3 + 4*b4] [ b1 + b2 - b3 + b4] [ b1 + b2 + b3 + 4*b5] [ -b2 + 2*b3 - b4 - b5 + b6] B.rescale_row(5, -1); B[0:6, 5] [-1/4*b1 - 5/4*b2 + 11/4*b3 - 2*b4] [ 3*b2 - 8*b3 + 6*b4] [-3/2*b1 + 3/2*b2 - 13/2*b3 + 4*b4] [ b1 + b2 - b3 + b4] [ b1 + b2 + b3 + 4*b5] [ b2 - 2*b3 + b4 + b5 - b6] B.add_multiple_of_row(4, 5,-1); B[0:6, 5] [-1/4*b1 - 5/4*b2 + 11/4*b3 - 2*b4] [ 3*b2 - 8*b3 + 6*b4] [-3/2*b1 + 3/2*b2 - 13/2*b3 + 4*b4] [ b1 + b2 - b3 + b4] [ b1 + 3*b3 - b4 + 3*b5 + b6] [ b2 - 2*b3 + b4 + b5 - b6] Notice that the last two entries of the final column now have just a single b1 and a single B2. We could continue to perform more row operations by hand, using the last two rows to progressively eliminate b1 and b2 from the other four expressions of the last column. Since the two last rows have zeros in their first five entries, only the entries in the final column would change. You will see that much of this section is about how to automate these final calculations. EEF Extended Echelon Form Sage will compute the extended echelon form, an improvement over what we did “by hand” in Sage MISLE. And the output can be requested as a “subdivided” matrix so that we can easily work with the pieces C and L. Pieces of subdivided matrices can be extracted with indices entirely similar to how we index the actual entries of a matrix. We will redo Example FS2 as an illustration of Theorem FS and Theorem PEEF. A = matrix(QQ,[[ 10, 0, 3, 8, 7], [-16, -1, -4, -10, -13], [ -6, 1, -3, -6, -6], [ 0, 2, -2, -3, -2], [ 3, 0, 1, 2, 3], [ -1, -1, 1, 1, 0]]) N = A.extended_echelon_form(subdivide=True); N [ 1 0 0 0 2| 0 0 1 -1 2 -1] [ 0 1 0 0 -3| 0 0 -2 3 -3 3] [ 0 0 1 0 1| 0 0 1 1 3 3] [ 0 0 0 1 -2| 0 0 -2 1 -4 0] Chapter M 70 [--------------+-----------------] [ 0 0 0 0 0| 1 0 3 -1 3 1] [ 0 0 0 0 0| 0 1 -2 1 1 -1] C = N.subdivision(0,0); C [ 1 0 0 0 2] [ 0 1 0 0 -3] [ 0 0 1 0 1] [ 0 0 0 1 -2] L = N.subdivision(1,1); L [ 1 0 3 -1 3 1] [ 0 1 -2 1 1 -1] K = N.subdivision(0,1); K [ 0 0 1 -1 2 -1] [ 0 0 -2 3 -3 3] [ 0 0 1 1 3 3] [ 0 0 -2 1 -4 0] C.rref() == C True L.rref() == L True A.right_kernel() == C.right_kernel() True A.row_space() == C.row_space() True A.column_space() == L.right_kernel() True Chapter M 71 A.left_kernel() == L.row_space() True J = N.matrix_from_columns(range(5, 11)); J [ 0 0 1 -1 2 -1] [ 0 0 -2 3 -3 3] [ 0 0 1 1 3 3] [ 0 0 -2 1 -4 0] [ 1 0 3 -1 3 1] [ 0 1 -2 1 1 -1] J.is_singular() False J*A [ 1 0 0 0 2] [ 0 1 0 0 -3] [ 0 0 1 0 1] [ 0 0 0 1 -2] [ 0 0 0 0 0] [ 0 0 0 0 0] J*A == A.rref() True Notice how we can employ the uniqueness of reduced row-echelon form (Theorem RREFU) to verify that C and L are in reduced row-echelon form. Realize too, that subdivided output is an optional behavior of the .extended_echelon_form() method and must be requested with subdivide=True. Additionally, it is the uniquely determined matrix J that provides us with a standard version of the final column of the row-reduced augmented matrix using variables in the vector of constants, as first shown back in Example CSANS and ver- ified here with variables defined as in Sage RRSM. (The vector method .column() converts a vector into a 1-column matrix, used here to format the matrix-vector product nicely.) R.<b1,b2,b3,b4,b5,b6> = QQ[] const = vector(R, [b1, b2, b3, b4, b5, b6]) (J*const).column() [ b3 - b4 + 2*b5 - b6] [-2*b3 + 3*b4 - 3*b5 + 3*b6] [ b3 + b4 + 3*b5 + 3*b6] [ -2*b3 + b4 - 4*b5] [b1 + 3*b3 - b4 + 3*b5 + b6] [ b2 - 2*b3 + b4 + b5 - b6] Chapter M 72 Chapter VS Vector Spaces Section S: Subspaces VS Vector Spaces Our conception of a vector space has become much broader with the introduction of abstract vector spaces — those whose elements (“vectors”) are not just column vectors, but polynomials, matrices, sequences, functions, etc. Sage is able to perform computations using many different abstract and advanced ideas (such as derivatives of functions), but in the case of linear algebra, Sage will primarily stay with vector spaces of column vectors. Chapter R, and specifically, Section VR and Sage SUTH2 will show us that this is not as much of a limitation as it might first appear. While limited to vector spaces of column vectors, Sage has an impressive range of capabilities for vector spaces, which we will detail throughout this chapter. You may have already noticed that many questions about abstract vector spaces can be translated into questions about column vectors. This theme will continue, and Sage commands we already know will often be helpful in answering these questions. Theorem SSS, Theorem NSMS, Theorem CSMS, Theorem RSMS and Theorem LNSMS each tells us that a certain set is a subspace. The first is the abstract version of creating a subspace via the span of a set of vectors, but still applies to column vectors as a special case. The remaining four all begin with a matrix and create a subspace of column vectors. We have created these spaces many times already, but notice now that the description Sage outputs explicitly says they are vector spaces, and that there are still some parts of the output that we need to explain. Here are two reminders, first a span, and then a vector space created from a matrix. V = QQ^4 v1 = vector(QQ, [ 1, -1, 2, 4]) v2 = vector(QQ, [-3, 0, 2, 1]) v3 = vector(QQ, [-1, -2, 6, 9]) W = V.span([v1, v2, v3]) W Vector space of degree 4 and dimension 2 over Rational Field Basis matrix: [ 1 0 -2/3 -1/3] [ 0 1 -8/3 -13/3] A = matrix(QQ, [[1, 2, -4, 0, -4], [0, 1, -1, -1, -1], [3, 2, -8, 4, -8]]) W = A.column_space() W 73 Chapter VS 74 Vector space of degree 3 and dimension 2 over Rational Field Basis matrix: [ 1 0 3] [ 0 1 -4] Section B: Bases B Bases Every vector space in Sage has a basis — you can obtain this with the vector space method .basis(), and the result is a list of vectors. Another method for a vector space is .basis_matrix() which outputs a matrix whose rows are the vectors of a basis. Sometimes one form is more convenient than the other, but notice that the description of a vector space chooses to print the basis matrix (since its display is just a bit easier to read). A vector space typically has many bases (infinitely many), so which one does Sage use? You will notice that the basis matrices displayed are in reduced row-echelon form — this is the defining property of the basis chosen by Sage. Here is Example RSB again as an example of how bases are provided in Sage. V = QQ^3 v1 = vector(QQ, [ 2, -3, 1]) v2 = vector(QQ, [ 1, 4, 1]) v3 = vector(QQ, [ 7, -5, 4]) v4 = vector(QQ, [-7, -6, -5]) W = V.span([v1, v2, v3, v4]) W Vector space of degree 3 and dimension 2 over Rational Field Basis matrix: [ 1 0 7/11] [ 0 1 1/11] W.basis() [ (1, 0, 7/11), (0, 1, 1/11) ] W.basis_matrix() [ 1 0 7/11] [ 0 1 1/11] SUTH0 Sage Under The Hood, Round 0 Or perhaps, “under the bonnet” if you learned your English in the Commonwealth. This is the first in a series that aims to explain how our knowledge of linear algebra theory helps us understand the design, construction and informed use of Sage. How does Sage determine if two vector spaces are equal? Especially since these are infinite sets? One approach would be to take a spanning set for the first vector space (maybe a minimal spanning set), and ask if each element of the spanning set is an element of the second vector space. If so, the first vector space is a subset of Chapter VS 75 the second. Then we could turn it around, and determine if the second vector space is a subset of the first. By Definition SE, the two vector spaces would be equal if both subset tests succeeded. However, each time we would test if an element of a spanning set lives in a second vector space, we would need to solve a linear system. So for two large vector spaces, this could take a noticeable amount of time. There is a better way, made possible by exploiting two important theorems. For every vector space, Sage creates a basis that uniquely identifies the vector space. We could call this a “canonical basis.” By Theorem REMRS we can span the row space of a matrix by the rows of any row-equivalent matrix. So if we begin with a vector space described by any basis (or any spanning set, for that matter), we can make a matrix with these rows as vectors, and the vector space is now the row space of the matrix. Of all the possible row-equivalent matrices, which would you pick? Of course, the reduced row-echelon version is useful, and here it is critical to realize this version is unique (Theorem RREFU). So for every vector space, Sage takes a spanning set, makes its vectors the rows of a matrix, row-reduces the matrix and tosses out the zero rows. The result is what Sage calls an “echelonized basis.” Now, two vector spaces are equal if, and only if, they have equal “echelonized basis matrices.” It takes some computation to form the echelonized basis, but once built, the comparison of two echelonized bases can proceed very quickly by perhaps just comparing entries of the echelonized basis matrices. You might create a vector space with a basis you prefer (a “user basis”), but Sage always has an echelonized basis at hand. If you do not specify some alter- nate basis, this is the basis Sage will create and provide for you. We can now continue a discussion we began back in Sage SSNS. We have consistently used the basis=’pivot’ keyword when we construct null spaces. This is because we initially prefer to see the basis described in Theorem BNS, rather than Sage’s default basis, the echelonized version. But the echelonized version is always present and available. A = matrix(QQ, [[14, -42, -2, -44, -42, 100, -18], [-40, 120, -6, 129, 135, -304, 28], [11, -33, 0, -35, -35, 81, -11], [-21, 63, -4, 68, 72, -161, 13], [-4, 12, -1, 13, 14, -31, 2]]) K = A.right_kernel(basis=’pivot’) K.basis_matrix() [ 3 1 0 0 0 0 0] [ 0 0 1 -1 1 0 0] [-1 0 -1 2 0 1 0] [ 1 0 -2 0 0 0 1] K.echelonized_basis_matrix() [ 1 0 0 0 -4 -2 -1] [ 0 1 0 0 12 6 3] [ 0 0 1 0 -2 -1 -1] [ 0 0 0 1 -3 -1 -1] NME5 Nonsingular Matrix Equivalences, Round 5 We can easily illustrate our latest equivalence for nonsingular matrices. A = matrix(QQ, [[ 2, 3, -3, 2, 8, -4], [ 3, 4, -4, 4, 8, 1], Chapter VS 76 [-2, -2, 3, -3, -2, -7], [ 0, 1, -1, 2, 3, 4], [ 2, 1, 0, 1, -4, 4], [ 1, 2, -2, 1, 7, -5]]) not A.is_singular() True V = QQ^6 cols = A.columns() V == V.span(cols) True V.linear_dependence(cols) == [] True C Coordinates For vector spaces of column vectors, Sage can quickly determine the coordinates of a vector relative to a basis, as guaranteed by Theorem VRRB. We illustrate some new Sage commands with a simple example and then apply them to orthonormal bases. The vectors v1 and v2 are linearly independent and thus span a subspace with a basis of size 2. We first create this subspace and let Sage determine the basis, then we illustrate a new vector space method, .subspace_with_basis(), that allows us to specify the basis. (This method is very similar to .span_of_basis(), except it preserves a subspace relationship with the original vector space.) Notice how the description of the vector space makes it clear that W has a user-specified basis. Notice too that the actual subspace created is the same in both cases. V = QQ^3 v1 = vector(QQ,[ 2, 1, 3]) v2 = vector(QQ,[-1, 1, 4]) U=V.span([v1,v2]) U Vector space of degree 3 and dimension 2 over Rational Field Basis matrix: [ 1 0 -1/3] [ 0 1 11/3] W = V.subspace_with_basis([v1, v2]) W Vector space of degree 3 and dimension 2 over Rational Field User basis matrix: [ 2 1 3] [-1 1 4] Chapter VS 77 U == W True Now we manufacture a third vector in the subspace, and request a coordinatiza- tion in each vector space, which has the effect of using a different basis in each case. The vector space method .coordinate_vector(v) computes a vector whose entries express v as a linear combination of basis vectors. Verify for yourself in each case below that the components of the vector returned really give a linear combination of the basis vectors that equals v3. v3 = 4*v1 + v2; v3 (7, 5, 16) U.coordinate_vector(v3) (7, 5) W.coordinate_vector(v3) (4, 1) Now we can construct a more complicated example using an orthonormal basis, specifically the one from Example CROB4, but we will compute over QQbar, the field of algebraic numbers. We form the four vectors of the orthonormal basis, install them as the basis of a vector space and then ask for the coordinates. Sage treats the square roots in the scalars as “symbolic” expressions, so we need to explicitly coerce them into QQbar before computing the scalar multiples. V = QQbar^4 x1 = vector(QQbar, [ 1+I, 1, 1-I, I]) x2 = vector(QQbar, [ 1+5*I, 6+5*I, -7-I, 1-6*I]) x3 = vector(QQbar, [-7+34*I, -8-23*I, -10+22*I, 30+13*I]) x4 = vector(QQbar, [ -2-4*I, 6+I, 4+3*I, 6-I]) v1 = QQbar(1/sqrt(6)) * x1 v2 = QQbar(1/sqrt(174)) * x2 v3 = QQbar(1/sqrt(3451))* x3 v4 = QQbar(1/sqrt(119)) * x4 W = V.subspace_with_basis([v1, v2, v3, v4]) w = vector(QQbar, [2, -3, 1, 4]) c = W.coordinate_vector(w); c (0.?e-14 - 2.0412414523194?*I, -1.4403862828000? + 2.27429413073671?*I, 2.0427196489446? - 3.59178204939423?*I, 0.55001909821693? + 1.10003819643385?*I) Is this right? Our exact coordinates in the text are printed differently, but we can check that they are the same numbers: c[0] == 1/sqrt(6)*(-5*I) Chapter VS 78 True c[1] == 1/sqrt(174)*(-19+30*I) True c[2] == 1/sqrt(3451)*(120-211*I) True c[3] == 1/sqrt(119)*(6+12*I) True With an orthonormal basis, we can illustrate Theorem CUMOS by making the four vectors the columns of 4 × 4 matrix and verifying the result is a unitary matrix. U = column_matrix([v1, v2, v3, v4]) U.is_unitary() True We will see coordinate vectors again, in a more formal setting, in Sage VR. Section D: Dimension D Dimension Now we recognize that every basis has the same size, even if there are many different bases for a given vector space. The dimension is an important piece of information about a vector space, so Sage routinely provides this as part of the description of a vector space. But it can be returned by itself with the vector space method .dimension(). Here is an example of a subspace with dimension 2. V = QQ^4 v1 = vector(QQ, [2, -1, 3, 1]) v2 = vector(QQ, [3, -3, 4, 0]) v3 = vector(QQ, [1, -2, 1, -1]) v4 = vector(QQ, [4, -5, 5, -1]) W = span([v1, v2, v3, v4]) W Vector space of degree 4 and dimension 2 over Rational Field Basis matrix: [ 1 0 5/3 1] [ 0 1 1/3 1] W.dimension() Chapter VS 79 2 RNM Rank and Nullity of a Matrix The rank and nullity of a matrix in Sage could be exactly what you would have guessed. But we need to be careful. The rank is the rank. But nullity in Sage is the dimension of the left nullspace. So we have matrix methods .nullity(), .left_nullity(), .right_nullity(), where the first two are equal and correspond to Sage’s preference for rows, and the third is the column version used by the text. That said, a “row version” of Theorem RPNC is also true. A = matrix(QQ, [[-1, 0, -4, -3, 1, -1, 0, 1, -1], [ 1, 1, 6, 6, 5, 3, 4, -5, 3], [ 2, 0, 7, 5, -3, 1, -1, -1, 2], [ 2, 1, 6, 6, 3, 1, 3, -3, 5], [-2, 0, -1, -1, 3, 3, 1, -3, -4], [-1, 1, 4, 4, 7, 5, 4, -7, -1]]) A.rank() 4 A.right_nullity() 5 A.rank() + A.right_nullity() == A.ncols() True A.rank() + A.left_nullity() == A.nrows() True NME6 Nonsingular Matrix Equivalences, Round 6 Recycling the nonsingular matrix from Sage NME5 we can use Sage to verify the two new equivalences of Theorem NME6. A = matrix(QQ, [[ 2, 3, -3, 2, 8, -4], [ 3, 4, -4, 4, 8, 1], [-2, -2, 3, -3, -2, -7], [ 0, 1, -1, 2, 3, 4], [ 2, 1, 0, 1, -4, 4], [ 1, 2, -2, 1, 7, -5]]) not A.is_singular() True A.rank() == A.ncols() Chapter VS 80 True A.right_nullity() == 0 True Section PD: Properties of Dimension DMS Dimensions of Matrix Subspaces The theorems in this section about the dimensions of various subspaces associated with matrices can be tested easily in Sage. For a large arbitrary matrix, we first verify Theorem RMRT, followed by the four conclusions of Theorem DFS. A = matrix(QQ, [ [ 1, -2, 3, 2, 0, 2, 2, -1, 3, 8, 0, 7], [-1, 2, -2, -1, 3, -3, 5, -2, -6, -7, 6, -2], [ 0, 0, 1, 1, 0, 0, 1, -3, -2, 0, 3, 8], [-1, -1, 0, -1, -1, 0, -6, -2, -5, -6, 5, 1], [ 1, -3, 2, 1, -4, 4, -6, 2, 7, 7, -5, 2], [-2, 2, -2, -2, 3, -3, 6, -1, -8, -8, 7, -7], [ 0, -3, 2, 0, -3, 3, -7, 1, 2, 3, -1, 0], [ 0, -1, 2, 1, 2, 0, 4, -3, -3, 2, 6, 6], [-1, 1, 0, -1, 2, -1, 6, -2, -6, -3, 8, 0], [ 0, -4, 4, 0, -2, 4, -4, -2, -2, 4, 8, 6] ]) m = A.nrows() n = A.ncols() r = A.rank() m, n, r (10, 12, 7) A.transpose().rank() == r True A.right_kernel().dimension() == n - r True A.column_space().dimension() == r True A.row_space().dimension() == r Chapter VS 81 True A.left_kernel().dimension() == m - r True Chapter D Determinants Section DM: Determinant of a Matrix EM Elementary Matrices Each of the three types of elementary matrices can be constructed easily, with syntax similar to methods for performing row operations on matrices. Here we have three 4 × 4 elementary matrices, in order: E1,3 , E2 (7), E2,4 (9). Notice the change in numbering on the rows, and the order of the parameters. A = elementary_matrix(QQ, 4, row1=0, row2=2); A [0 0 1 0] [0 1 0 0] [1 0 0 0] [0 0 0 1] A = elementary_matrix(QQ, 4, row1=1, scale=7); A [1 0 0 0] [0 7 0 0] [0 0 1 0] [0 0 0 1] A = elementary_matrix(QQ, 4, row1=3, row2=1, scale=9); A [1 0 0 0] [0 1 0 0] [0 0 1 0] [0 9 0 1] Notice that row1 is always the row that is being changed. The keywords can be removed, but the scale keyword must be used to create the second type of elementary matrix, to avoid confusion with the first type. We can illustrate some of the results of this section with two examples. First, we convert a matrix into a second matrix which is row-equivalent, and then accomplish the same thing with matrix multiplication and a product of elementary matrices. A = matrix(QQ, [[6, -2, 3, -2], [3, 3, 1, 8], [4, 0, 5, 4]]) B = copy(A) B.swap_rows(0,2) 82 Chapter D 83 E1 = elementary_matrix(QQ, 3, row1=0, row2=2) B.rescale_row(1, 5) E2 = elementary_matrix(QQ, 3, row1=1, scale=5) B.add_multiple_of_row(1, 0, -3) E3 = elementary_matrix(QQ, 3, row1=1, row2=0, scale=-3) B [ 4 0 5 4] [ 3 15 -10 28] [ 6 -2 3 -2] E3*E2*E1*A [ 4 0 5 4] [ 3 15 -10 28] [ 6 -2 3 -2] B == E3*E2*E1*A True R = E3*E2*E1 R.is_singular() False R [ 0 0 1] [ 0 5 -3] [ 1 0 0] The matrix R, the product of three elementary matrices, can be construed as the collective effect of the three row operations employed. With more row operations, R might look even less like an identity matrix. As the product of nonsingular matrices (Theorem EMN), R is nonsingular (Theorem NPNT). The matrix B above is not in reduced row-echelon form (it was just row-equivalent to A). What if we were to begin with a matrix and track all of the row operations required to bring the matrix to reduced row-echelon form? As above, we could form the associated elementary matrices and form their product, creating a single matrix R that captures all of the row operations. It turns out we have already done this. Extended echelon form is the subject of Theorem PEEF, whose second conclusion says that B = JA, where A is the original matrix, and B is the row-equivalent matrix in reduced row-echelon form. Then J is a square nonsingular matrix that is the product of the sequence of elementary matrices associated with the sequence of row operations converting A into B. There may be many, many different sequences of row operations that convert A to B, but the requirement that extended echelon form be in reduced row-echelon form guarantees that J is unique. DM Determinant of a Matrix Computing the determinant in Sage is straightforward. Chapter D 84 A = matrix(QQ,[[ 4, 3, -2, -9, -11, -14, -4, 11, -4, 4], [ 0, 1, 0, -1, -1, -3, -2, 3, 6, 15], [ 0, 1, 1, 1, 3, 0, -3, -5, -9, -3], [-3, 3, 3, 6, 12, 3, -9, -12, -9, 15], [-2, 0, 2, 3, 5, 3, -3, -5, -3, 7], [ 3, 3, -1, -7, -8, -12, -5, 8, -4, 9], [ 0, -1, 1, 2, 2, 5, 3, -6, -8, -12], [-1, -2, -1, 1, -2, 5, 6, 3, 6, -8], [ 2, 1, -3, -6, -8, -7, 0, 12, 11, 8], [-3, -2, 2, 5, 5, 9, 2, -7, -4, -3]]) A.determinant() 3 Random matrices, even with small entries, can have very large determinants. B = matrix(QQ, 15, [ [-5, -5, -5, 4, 1, 1, -3, 0, 4, 4, -2, -4, 2, 3, -1], [ 1, 1, -4, 3, 3, 4, 1, -1, -5, 4, -3, 0, -1, 0, 0], [ 3, 4, -2, 3, -1, -5, -1, -4, -5, 0, -1, 2, -4, -1, -1], [ 2, -4, 4, -3, -3, -3, -1, -3, -3, -1, 2, 4, -1, -1, -3], [ 1, 4, -3, -1, -2, -2, 1, -1, 3, -5, -4, -2, -2, -2, -5], [-1, -2, 4, 0, 4, 1, 1, 4, -5, 3, 1, -1, 4, 2, -2], [ 4, 3, 2, 4, 4, -5, 2, -5, -5, 2, -5, -4, -4, 0, 3], [ 1, -2, 0, -2, -2, 0, 2, 3, 1, 2, -4, 0, -5, -2, 2], [-4, -4, 2, 1, -1, 4, -2, 1, -2, 2, -1, -1, 3, 4, -1], [-4, 0, 2, 3, -4, -5, 3, -5, 4, -4, -2, 3, 3, -3, 0], [ 1, 2, 3, -4, 2, 0, -4, -1, 1, -3, 1, 4, -2, 4, -1], [-3, 3, 0, 2, 1, -2, -4, 0, -1, -1, -1, 2, 3, 1, -4], [ 4, 3, -3, -4, 3, 1, -3, 2, 1, -5, -5, -3, 2, 1, 4], [-3, -5, -1, -5, -2, 0, -3, 1, 2, -1, 0, -4, 3, -2, 3], [-1, 1, -3, -1, 3, -3, 2, -3, -5, -1, -1, 3, -1, 2, 3] ]) B.determinant() 202905135564 Sage is incredibly fast at computing determinants with rational entries. Try the following two compute cells on whatever computer you might be using right now. The one unfamilar command clears the value of the determinant that Sage caches, so we get accurate timing information across multiple evaluations. C = random_matrix(QQ, 100, 100) timeit("C.determinant(); C._cache={}") # random 5 loops, best of 3: 152 ms per loop C.determinant() # random -54987836999...175801344 Chapter D 85 Section PDM: Properties of Determinants of Ma- trices NME7 Nonsingular Matrix Equivalences, Round 7 Our latest addition to the NME series is sometimes considered the best compu- tational way to determine if a matrix is nonsingular. While computing a single number for a comparison to zero is a very appealing process, other properties of nonsingular matrices are faster for determining if a matrix is nonsingular. Still, it is a worthwhile equivalence. We illustrate both directions. A = matrix(QQ, [[ 1, 0, 2, 1, 1, -3, 4, -6], [ 0, 1, 0, 5, 3, -8, 0, 6], [ 1, 1, 4, 4, 2, -8, 4, -2], [ 0, 1, 6, 0, -2, -1, 1, 0], [-1, 1, 0, 1, 0, -2, -1, 5], [ 0, -1, -6, -3, 1, 4, 1, -5], [-1, -1, -2, -3, -2, 7, -4, 2], [ 0, 2, 0, 0, -2, -2, 0, 6]]) A.is_singular() False A.determinant() 4 B = matrix(QQ, [[ 2, -1, 0, 4, -6, -2, 8, -2], [-1, 1, -1, -4, -1, 7, -7, -2], [-1, 1, 1, 0, 7, -5, -6, 7], [ 2, -1, -1, 3, -5, -2, 5, 2], [ 1, -2, 0, 2, -1, -3, 8, 0], [-1, 1, 0, -2, 6, -1, -8, 5], [ 2, -1, -1, 2, -7, 1, 7, -3], [-1, 0, 1, 0, 6, -4, -2, 4]]) B.is_singular() True B.determinant() 0 PDM Properties of Determinants of Matrices There are many properties of determinants detailed in this section. You have the tools in Sage to explore all of them. We will illustrate with just two. One is not par- ticularly deep, but we personally find the interplay between matrix multiplication and the determinant nothing short of amazing, so we will not pass up the chance to verify Theorem DRMM. Add some checks of other properties yourself. We copy the sixth row of the matrix C to the fourth row, so by Theorem DERC, the determinant is zero. Chapter D 86 C = matrix(QQ, [[-3, 4, 0, 1, 2, 0, 5, 0], [ 4, 0, -4, 7, 0, 5, 0, 5], [ 7, 4, 2, 2, 4, 0, -2, 8], [ 5, -4, 8, 2, 6, -1, -4, -4], [ 8, 0, 7, 4, 7, 5, 2, -3], [ 6, 5, 3, 7, 4, 2, 4, -3], [ 1, 6, -4, -4, 3, 8, 5, -2], [ 2, 4, -2, 8, 2, 5, 2, 2]]) C[3] = C[5] C [-3 4 0 1 2 0 5 0] [ 4 0 -4 7 0 5 0 5] [ 7 4 2 2 4 0 -2 8] [ 6 5 3 7 4 2 4 -3] [ 8 0 7 4 7 5 2 -3] [ 6 5 3 7 4 2 4 -3] [ 1 6 -4 -4 3 8 5 -2] [ 2 4 -2 8 2 5 2 2] C.determinant() 0 Now, a demonstration of Theorem DRMM. A = matrix(QQ, [[-4, 7, -2, 6, 8, 0, -4, 6], [ 1, 6, 5, 8, 2, 3, 1, -1], [ 5, 0, 7, -3, 7, -3, 6, -3], [-4, 5, 8, 3, 6, 8, -1, -1], [ 0, 0, -3, -3, 4, 4, 2, 5], [-3, -2, -3, 8, 8, -3, -1, 1], [-4, 0, 2, 4, 4, 4, 7, 2], [ 3, 3, -4, 5, -2, 3, -1, 5]]) B = matrix(QQ, [[ 1, 2, -4, 8, -1, 2, -1, -3], [ 3, 3, -2, -4, -3, 8, 1, 6], [ 1, 8, 4, 0, 4, -2, 0, 8], [ 6, 8, 1, -1, -4, -3, -2, 5], [ 0, 5, 1, 4, -3, 2, -3, -2], [ 2, 4, 0, 7, 8, -1, 5, 8], [ 7, 1, 1, -1, -1, 7, -2, 6], [ 2, 3, 4, 7, 3, 4, 7, -2]]) C = A*B; C [ 35 99 28 12 -51 46 25 16] [ 83 144 11 -3 -17 20 -11 141] [ 24 62 6 23 -25 52 -68 30] [ 44 153 42 19 53 0 20 169] [ 11 5 11 80 33 53 45 -13] [ 25 58 23 -5 -79 -24 -45 -35] [ 83 89 47 15 15 37 4 110] [ 47 39 -12 56 -2 29 48 14] Chapter D 87 Adet = A.determinant(); Adet 9350730 Bdet = B.determinant(); Bdet 6516260 Cdet = C.determinant(); Cdet 60931787869800 Adet*Bdet == Cdet True Earlier, in Sage NME1 we used the random_matrix() constructor with the algorithm=’unimodular’ keyword to create a nonsingular matrix. We can now reveal that in this context, “unimodular” means “with determinant 1.” So such a matrix will always be nonsingular by Theorem SMZD. But more can be said. It is not obvious at all, and Solution PDM.M20 has just a partial explanation, but the inverse of a unimodular matrix with all integer entries will have an inverse with all integer entries. With a fraction-free inverse many “textbook” exercises can be constructed through the use of unimodular matrices. Experiment for yourself below. An upper_bound keyword can be used to control the size of the entries of the matrix, though this will have little control over the inverse. A = random_matrix(QQ, 5, algorithm=’unimodular’) A, A.inverse() # random ( [ -9 -32 -118 273 78] [-186 30 22 -18 375] [ 2 9 31 -69 -21] [ 52 -8 -8 3 -105] [ 4 15 54 -120 -39] [-147 25 19 -12 297] [ -3 -11 -38 84 28] [ -47 8 6 -4 95] [ -5 -18 -66 152 44], [ -58 10 7 -5 117] ) Chapter E Eigenvalues Section EE: Eigenvalues and Eigenvectors EE Eigenvalues and Eigenvectors Sage can compute eigenvalues and eigenvectors of matrices. We will see shortly that there are subtleties involved with using these routines, but here is a quick example to begin with. These two commands should be enough to get you started with most of the early examples in this section. See the end of the section for more comprehensive advice. For a square matrix, the methods .eigenvalues() and .eigenvectors_right() will produce what you expect, though the format of the eigenvector output requires some explanation. Here is Example SEE from the start of this chapter. A = matrix(QQ, [[ 204, 98, -26, -10], [-280, -134, 36, 14], [ 716, 348, -90, -36], [-472, -232, 60, 28]]) A.eigenvalues() [4, 0, 2, 2] A.eigenvectors_right() [(4, [ (1, -1, 2, 5) ], 1), (0, [ (1, -4/3, 10/3, -4/3) ], 1), (2, [ (1, 0, 7, 2), (0, 1, 3, 2) ], 2)] The three eigenvalues we know are included in the output of eigenvalues(), though for some reason the eigenvalue λ = 2 shows up twice. The output of the eigenvectors_right() method is a list of triples. Each triple begins with an eigenvalue. This is followed by a list of eigenvectors for that eigen- value. Notice the first eigenvector is identical to the one we described in Example SEE. The eigenvector for λ = 0 is different, but is just a scalar multiple of the one from Example SEE. For λ = 2, we now get two eigenvectors, and neither looks like either of the ones from Example SEE. (Hint: try writing the eigenvectors from the example as linear combinations of the two in the Sage output.) An explanation of the the third part of each triple (an integer) will have to wait, though it can be optionally surpressed if desired. 88 Chapter E 89 One cautionary note: The word lambda has a special purpose in Sage, so do not try to use this as a name for your eigenvalues. CEVAL Computing Eigenvalues We can now give a more careful explantion about eigenvalues in Sage. Sage will com- pute the characteristic polynomial of a matrix, with amazing ease (in other words, quite quickly, even for large matrices). The two matrix methods .charpoly() and .characteristic_polynomial() do exactly the same thing. We will use the longer name just to be more readable, you may prefer the shorter. We now can appreciate a very fundamental obstacle to determining the eigen- values of a matrix, which is a theme that will run through any advanced study of linear algebra. Study this example carefully before reading the discussion that follows. A = matrix(QQ, [[-24, 6, 0, -1, 31, 7], [ -9, -2, -8, -17, 24, -29], [ 4, -10, 1, 1, -12, -36], [-19, 11, -1, -4, 33, 29], [-11, 6, 2, 3, 14, 21], [ 5, -1, 2, 5, -11, 4]]) A.characteristic_polynomial() x^6 + 11*x^5 + 15*x^4 - 84*x^3 - 157*x^2 + 124*x + 246 A.characteristic_polynomial().factor() (x + 3) * (x^2 - 2) * (x^3 + 8*x^2 - 7*x - 41) B = A.change_ring(QQbar) B.characteristic_polynomial() x^6 + 11*x^5 + 15*x^4 - 84*x^3 - 157*x^2 + 124*x + 246 B.characteristic_polynomial().factor() (x - 2.356181157637288?) * (x - 1.414213562373095?) * (x + 1.414213562373095?) * (x + 2.110260216209409?) * (x + 3) * (x + 8.24592094142788?) We know by Theorem EMRCP that to compute eigenvalues, we need the roots of the characteristic polynomial, and from basic algebra, we know these correspond to linear factors. However, with our matrix defined with entries from QQ, the factor- ization of the characteristic polynomial does not “leave” that number system, only factoring “far enough” to retain factors √ with rational coefficients. The solutions to x2 − 2 = 0 are somewhat obvious (± 2 ≈ ±1.414213), but the roots of the cubic factor are more obscure. But then we have QQbar to the rescue. Since this number system contains the roots of every possible polynomial with integer coefficients, we can totally factor any characteristic polynomial that results from a matrix with entries from QQbar. A common situation will be to begin with a matrix having rational entries, yet the matrix has a characteristic polynomial with roots that are complex numbers. We can demonstrate this behavior with the extend keyword option, which tells Sage whether or not to expand the number system to contain the eigenvalues. Chapter E 90 A.eigenvalues(extend=False) [-3] A.eigenvalues(extend=True) [-3, -1.414213562373095?, 1.414213562373095?, -8.24592094142788?, -2.110260216209409?, 2.356181157637288?] For matrices with entries from QQ, the default behavior is to extend to QQbar when necessary. But for other number systems, you may need to explicitly use the extend=True option. From a factorization of the characteristic polynomial, we can see the algebraic multiplicity of each eigenvalue as the second entry of the each pair returned in the list. We demonstrate with Example SEE, extending to QQbar, which is not strictly necessary for this simple matrix. A = matrix(QQ, [[204, 98, -26, -10], [-280, -134, 36, 14], [716, 348, -90, -36], [-472, -232, 60, 28]]) A.characteristic_polynomial().roots(QQbar) [(0, 1), (2, 2), (4, 1)] One more example, which illustrates the behavior when we use floating-point approximations as entries (in other words, we use CDF as our number system). This is Example EMMS4, both as an exact matrix with entries from QQbar and as an approximate matrix with entries from CDF. A = matrix(QQ, [[-2, 1, -2, -4], [12, 1, 4, 9], [ 6, 5, -2, -4], [ 3, -4, 5, 10]]) A.eigenvalues() [1, 2, 2, 2] B = A.change_ring(CDF) B.eigenvalues() [2.00001113832 + 1.75245331504e-05*I, 2.00000960797 - 1.84081096655e-05*I, 1.99997925371 + 8.8357651007e-07*I, 1.0] So, we see λ = 2 as an eigenvalue with algebraic multiplicity 3, while the nu- merical results contain three complex numbers, each very, very close to 2. The ap- proximate nature of these eigenvalues may be disturbing (or alarming). However, their computation, as floating-point numbers, can be incredibly fast with sophis- ticated algorithms allowing the analysis of huge matrices with millions of entries. And perhaps your original matrix includes data from an experiment, and is not even exact in the first place. Designing and analyzing algorithms to perform these Chapter E 91 computations quickly and accurately is part of the field known as numerical linear algebra. One cautionary note: Sage uses a definition of the characteristic polynomial slightly different than ours, namely det (xIn − A). This has the advantage that the xn term always has a positive one as the leading coefficient. For even values of n the two definitions create the identical polynomial, and for odd values of n, the two polynomials differ only by a multiple of −1. The reason this is not very critical is that Theorem EMRCP is true in either case, and this is a principal use of the characteristic polynomial. Our definition is more amenable to computations by hand. CEVEC Computing Eigenvectors There are three ways to get eigenvectors in Sage. For each eigenvalue, the method .eigenvectors_right() will return a list of eigenvectors that is a basis for the as- sociated eigenspace. The method .eigenspaces_right() will return an eigenspace (in other words, a vector space, rather than a list of vectors) for each eigenvalue. There are also eigenmatrix methods which we will describe at the end of the chapter in Sage MD. The matrix method .eigenvectors_right() (or equivalently the matrix method .right_eigenvectors()) produces a list of triples, one triple per eigenvalue. Each triple has an eigenvalue, a list, and then the algebraic multiplicity of the eigenvalue. The list contains vectors forming a basis for the eigenspace. Notice that the length of the list of eigenvectors will be the geometric multiplicity (and there is no easier way to get this information). A = matrix(QQ, [[-5, 1, 7, -15], [-2, 0, 5, -7], [ 1, -5, 7, -3], [ 4, -7, 3, 4]]) A.eigenvectors_right() [(3, [ (1, 0, -1, -1), (0, 1, 2, 1) ], 2), (-2*I, [(1, 0.3513513513513514? + 0.10810810810810811?*I, 0.02702702702702703? + 0.1621621621621622?*I, -0.2972972972972973? + 0.2162162162162163?*I)], 1), (2*I, [(1, 0.3513513513513514? - 0.10810810810810811?*I, 0.02702702702702703? - 0.1621621621621622?*I, -0.2972972972972973? - 0.2162162162162163?*I)], 1)] Note that this is a good place to practice burrowing down into Sage output that is full of lists (and lists of lists). See if you can extract just the second eigenvector for λ = 3 using a single statement. Or perhaps try obtaining the geometric multiplicity of λ = −2i. Notice, too, that Sage has automatically upgraded to QQbar to get the complex eigenvalues. The matrix method .eigenspaces_right() (equal to .right_eigenspaces()) produces a list of pairs, one pair per eigenvalue. Each pair has an eigenvalue, fol- lowed by the eigenvalue’s eigenspace. Notice that the basis matrix of the eigenspace may not have the same eigenvectors you might get from other methods. Similar to the eigenvectors method, the dimension of the eigenspace will yield the geometric multiplicity (and there is no easier way to get this information). If you need the alge- braic multiplicities, you can supply the keyword option algebraic_multiplicity=True to get back triples with the algebraic multiplicity in the third entry of the triple. We will recycle the example above, and not demonstrate the algebraic multiplicity option. (We have formatted the one-row basis matrices over QQbar across several lines.) Chapter E 92 A.eigenspaces_right() [ (3, Vector space of degree 4 and dimension 2 over Rational Field User basis matrix: [ 1 0 -1 -1] [ 0 1 2 1]), (-2*I, Vector space of degree 4 and dimension 1 over Algebraic Field User basis matrix: [ 1 0.3513513513513514? + 0.10810810810810811?*I 0.02702702702702703? + 0.1621621621621622?*I -0.2972972972972973? + 0.2162162162162163?*I]), (2*I, Vector space of degree 4 and dimension 1 over Algebraic Field User basis matrix: [ 1 0.3513513513513514? - 0.10810810810810811?*I 0.02702702702702703? - 0.1621621621621622?*I -0.2972972972972973? - 0.2162162162162163?*I]) ] Notice how the output includes a subspace of dimension two over the rationals, and two subspaces of dimension one over the algebraic numbers. The upcoming Subsection EE.ECEE has many examples, which mostly reflect techniques that might be possible to verify by hand. Here is the same matrix as above, analyzed in a similar way. Practicing the examples in this subsection, either directly with the higher-level Sage commands, and/or with more primitive commands (as below) would be an extremely good exercise at this point. A = matrix(QQ, [[-5, 1, 7, -15], [-2, 0, 5, -7], [ 1, -5, 7, -3], [ 4, -7, 3, 4]]) # eigenvalues A.characteristic_polynomial() x^4 - 6*x^3 + 13*x^2 - 24*x + 36 A.characteristic_polynomial().factor() (x - 3)^2 * (x^2 + 4) A.characteristic_polynomial().roots(QQbar) [(3, 2), (-2*I, 1), (2*I, 1)] # eigenspaces, rational (A-3*identity_matrix(4)).right_kernel(basis=’pivot’) Vector space of degree 4 and dimension 2 over Rational Field User basis matrix: [ 1 1 1 0] Chapter E 93 [-2 -1 0 1] # eigenspaces, complex B = A.change_ring(QQbar) (B-(2*I)*identity_matrix(4)).right_kernel(basis=’pivot’) Vector space of degree 4 and dimension 1 over Algebraic Field User basis matrix: [8/5*I - 11/5 4/5*I - 3/5 2/5*I + 1/5 1] (B-(-2*I)*identity_matrix(4)).right_kernel(basis=’pivot’) Vector space of degree 4 and dimension 1 over Algebraic Field User basis matrix: [-8/5*I - 11/5 -4/5*I - 3/5 -2/5*I + 1/5 1] Notice how we changed the number system to the algebraic numbers before working with the complex eigenvalues. Also, we are using the basis=’pivot’ keyword option so that bases for the eigenspaces look more like the bases described in Theorem BNS. By now, it should be clear why we keep using the “right” variants of these methods. Eigenvectors can be defined “on the right”, Ax = λx as we have done, or “on the left,” xA = λx. So use the “right” versions of the eigenvalue and eigenvector commands to stay consistent with the text. Recognize, too, that eigenspaces may be computed with different bases than those given in the text (typically like those for null spaces with the basis=’echelon’ option). Why does the .eigenvalues() method not come in left and right versions? The upcoming Theorem ETM can be used to show that the two versions would have identical output, so there is no need. Section PEE: Properties of Eigenvalues and Eigen- vectors NME8 Nonsingular Matrix Equivalences, Round 8 Zero eigenvalues are another marker of singular matrices. We illustrate with two matrices, the first nonsingular, the second not. A = matrix(QQ, [[ 1, 0, 2, 1, 1, -3, 4, -6], [ 0, 1, 0, 5, 3, -8, 0, 6], [ 1, 1, 4, 4, 2, -8, 4, -2], [ 0, 1, 6, 0, -2, -1, 1, 0], [-1, 1, 0, 1, 0, -2, -1, 5], [ 0, -1, -6, -3, 1, 4, 1, -5], [-1, -1, -2, -3, -2, 7, -4, 2], [ 0, 2, 0, 0, -2, -2, 0, 6]]) A.is_singular() False 0 in A.eigenvalues() Chapter E 94 False B = matrix(QQ, [[ 2, -1, 0, 4, -6, -2, 8, -2], [-1, 1, -1, -4, -1, 7, -7, -2], [-1, 1, 1, 0, 7, -5, -6, 7], [ 2, -1, -1, 3, -5, -2, 5, 2], [ 1, -2, 0, 2, -1, -3, 8, 0], [-1, 1, 0, -2, 6, -1, -8, 5], [ 2, -1, -1, 2, -7, 1, 7, -3], [-1, 0, 1, 0, 6, -4, -2, 4]]) B.is_singular() True 0 in B.eigenvalues() True Section SD: Similarity and Diagonalization SM Similar Matrices It is quite easy to determine if two matrices are similar, using the matrix method .is_similar(). However, computationally this can be a very difficult proposition, so support in Sage is incomplete now, though it will always return a result for matrices with rational entries. Here are examples where the two matrices are and are not similar. Notice that the keyword option transformation=True will cause a pair to be returned, such that if the matrices are indeed similar, the matrix effecting the similarity transformation will be in the second slot of the pair. A = matrix(QQ, [[ 25, 2, -8, -1, 11, 26, 35], [ 28, 2, -15, 2, 6, 34, 31], [ 1, -17, -25, 28, -44, 26, -23], [ 36, -2, -24, 10, -1, 50, 39], [ 0, -7, -13, 14, -21, 14, -11], [-22, -17, -16, 27, -51, 1, -53], [ -1, 10, 17, -18, 28, -18, 15]]) B = matrix(QQ, [[-89, -16, -55, -23, -104, -48, -67], [-15, 1, -20, -21, -20, -60, -26], [ 48, 6, 37, 25, 59, 64, 46], [-96, -16, -49, -16, -114, -23, -67], [ 56, 10, 33, 13, 67, 29, 37], [ 10, 2, 2, -2, 12, -9, 4], [ 28, 6, 13, 1, 32, -4, 16]]) is_sim, trans = A.is_similar(B, transformation=True) is_sim True Since we knew in advance these two matrices are similar, we requested the transformation matrix, so the output is a pair. The similarity matrix is a bit of a mess, so we will use three Sage routines to clean up trans. We convert the entries to numerical approximations, clip very small values (less than 10−5 ) to zero and then round to three decimal places. You can experiment printing just trans all by itself. Chapter E 95 trans.change_ring(RDF).zero_at(10^-5).round(3) [ 1.0 0.0 0.0 0.0 0.0 0.0 0.0] [ 0.0 1.188 0.375 -0.375 0.562 -0.375 0.375] [-3.107 -1.072 0.764 1.043 -2.288 -1.758 -2.495] [ 7.596 0.3 -1.026 -7.29 11.306 -2.705 -1.6] [ 0.266 0.079 -0.331 -0.387 0.756 0.447 0.429] [-2.157 0.071 0.469 1.893 -3.014 0.355 0.196] [-0.625 0.291 -0.068 0.995 -1.315 1.125 1.179] The matrix C is not similar to A (and hence not similar to B by Theorem SER), so we illustrate the return value when we do not request the similarity matrix (since it does not even exist). C = matrix(QQ, [[ 16, -26, 19, -56, 26, 8, -49], [ 20, -43, 36, -102, 52, 23, -65], [-18, 29, -21, 62, -30, -9, 56], [-17, 31, -27, 73, -37, -16, 49], [ 18, -36, 30, -82, 43, 18, -54], [-32, 62, -56, 146, -77, -35, 88], [ 11, -19, 17, -44, 23, 10, -29]]) C.is_similar(A) False MD Matrix Diagonalization The third way to get eigenvectors is the matrix method .eigenmatrix_right() (and the analogous .eigenmatrix_left()). It always returns two square matrices of the same size as the original matrix. The first matrix of the output is a diagonal matrix with the eigenvalues of the matrix filling the diagonal entries of the matrix. The second matrix has eigenvectors in the columns, in the same order as the cor- responding eigenvalues. For a single eigenvalue, these columns/eigenvectors form a linearly independent set. A careful reading of the previous paragraph suggests the question: what if we do not have enough eigenvectors to fill the columns of the second square matrix? When the geometric multiplicity does not equal the algebraic multiplicity, the deficit is met by inserting zero columns in the matrix of eigenvectors. Conversely, when the matrix is diagonalizable, by Theorem DMFE the geometric and algebraic multiplicities of each eigenvalue are equal, and the union of the bases of the eigenspaces provides a complete set of linearly independent vectors. So for a matrix A, Sage will output two matrices, D and S such that S −1 AS = D. We can rewrite the relation above as AS = SD. In the case of a non-diagonalizable matrix, the matrix of eigenvectors is singular (it has zero columns), but the rela- tionship AS = SD still holds. Here are examples of the two scenarios, along with demonstrations of the matrix method is_diagonalizable(). A = matrix(QQ, [[ 2, -18, -68, 64, -99, -87, 83], [ 4, -10, -41, 34, -58, -57, 46], [ 4, 16, 59, -60, 86, 70, -72], [ 2, -15, -65, 57, -92, -81, 78], [-4, -7, -32, 31, -45, -31, 41], [ 2, -6, -22, 20, -32, -31, 26], [ 0, 7, 30, -27, 42, 37, -36]]) D, S = A.eigenmatrix_right() D Chapter E 96 [ 0 0 0 0 0 0 0] [ 0 2 0 0 0 0 0] [ 0 0 2 0 0 0 0] [ 0 0 0 -1 0 0 0] [ 0 0 0 0 -1 0 0] [ 0 0 0 0 0 -3 0] [ 0 0 0 0 0 0 -3] S [ 1 1 0 1 0 1 0] [ 13/23 0 1 0 1 0 1] [-20/23 1/3 -2 -1 2/7 -4/3 5/6] [ 21/23 1/3 1 0 10/7 1 0] [ 10/23 -1/6 1 -1 15/7 4/3 -4/3] [ 8/23 1/3 0 1 -1 0 1/2] [-10/23 1/6 -1 -1 6/7 -1/3 -1/6] S.inverse()*A*S == D True A.is_diagonalizable() True Now for a matrix that is far from diagonalizable. B = matrix(QQ, [[ 37, -13, 30, 81, -74, -13, 18], [ 6, 26, 21, -11, -46, -48, 19], [ 16, 10, 29, 16, -42, -39, 26], [-24, 8, -24, -53, 54, 13, -15], [ -8, 3, -8, -20, 24, 4, -5], [ 31, 12, 46, 48, -97, -56, 35], [ 8, 5, 16, 12, -34, -20, 11]]) D, S = B.eigenmatrix_right() D [-2 0 0 0 0 0 0] [ 0 -2 0 0 0 0 0] [ 0 0 -2 0 0 0 0] [ 0 0 0 6 0 0 0] [ 0 0 0 0 6 0 0] [ 0 0 0 0 0 6 0] [ 0 0 0 0 0 0 6] S [ 1 0 0 1 0 0 0] [ 1/4 0 0 1/4 0 0 0] [ 1/2 0 0 3/8 0 0 0] [-3/4 0 0 -5/8 0 0 0] [-1/4 0 0 -1/4 0 0 0] [ 1 0 0 7/8 0 0 0] Chapter E 97 [ 1/4 0 0 1/4 0 0 0] B*S == S*D True B.is_diagonalizable() False Chapter LT Linear Transformations Section LT: Linear Transformations LTS Linear Transformations, Symbolic There are several ways to create a linear transformation in Sage, and many natural operations on these objects are possible. As previously, rather than use the complex numbers as our number system, we will instead use the rational numbers, since Sage can model them exactly. We will also use the following transformation repeatedly for examples, when no special properties are required: T : C3 → C4 −x1 + 2x3 x1 x1 + 3x2 + 7x3 T x2 = x1 + x2 + x3 x3 2x1 + 3x2 + 5x3 To create this linear transformation in Sage, we first create a “symbolic” func- tion, which requires that we also first define some symbolic variables which are x1, x2 and x3 in this case. (You can bypass the intermediate variable outputs in your own work. We will use it consistently to allow us to spread the definition across several lines without the Sage preparser getting in the way. In other words, it is safe to combine the the two lines below and not use outputs.) x1, x2, x3 = var(’x1, x2, x3’) outputs = [ -x1 + 2*x3, x1 + 3*x2 + 7*x3, x1 + x2 + x3, 2*x1 + 3*x2 + 5*x3] T_symbolic(x1, x2, x3) = outputs You can experiment with T_symbolic, evaluating it at triples of rational num- bers, and perhaps doing things like calculating its partial derivatives. We will use it as input to the linear_transformation() constructor. We just need to specify carefully the domain and codomain, now as vector spaces over the rationals rather than the complex numbers. T = linear_transformation(QQ^3, QQ^4, T_symbolic) You can now, of course, experiment with T via tab-completion, but we will explain the various properties of Sage linear transformations as we work through this chapter. Even some seemingly simple operations, such as printing T will require some explanation. But for starters, we can evaluate T. 98 Chapter LT 99 u = vector(QQ, [3, -1, 2]) w = T(u) w (1, 14, 4, 13) w.parent() Vector space of dimension 4 over Rational Field Here is a small verification of Theorem LTTZZ. zero_in = zero_vector(QQ, 3) zero_out = zero_vector(QQ, 4) T(zero_in) == zero_out True Note that Sage will recognize some symbolic functions as not being linear trans- formations (in other words, inconsistent with Definition LT), but this detection is fairly easy to fool. We will see some safer ways to create a linear transformation shortly. LTM Linear Transformations, Matrices A second way to build a linear transformation is to use a matrix, as motivated by Theorem MBLT. But there is one caveat. We have seen that Sage has a preference for rows, so when defining a linear transformation with a product of a matrix and a vector, Sage forms a linear combination of the rows of the matrix with the scalars of the vector. This is expressed by writing the vector on the left of the matrix, where if we were to interpret the vector as a 1-row matrix, then the definition of matrix multiplication would do the right thing. Remember throughout, a linear transformation has very little to do with the mechanism by which we define it. Whether or not we use matrix multiplication with vectors on the left (Sage internally) or matrix multiplication with vectors on the right (your text), what matters is the function that results. One concession to the “vector on the right” approach is that we can tell Sage that we mean for the matrix to define the linear transformation by multiplication with the vector on the right. Here is our running example again — with some explanation following. A = matrix(QQ, [[-1, 0, 2], [ 1, 3, 7], [ 1, 1, 1], [ 2, 3, 5]]) T = linear_transformation(QQ^3, QQ^4, A, side=’right’) T Vector space morphism represented by the matrix: [-1 1 1 2] [ 0 3 1 3] [ 2 7 1 5] Domain: Vector space of dimension 3 over Rational Field Codomain: Vector space of dimension 4 over Rational Field The way T prints reflects the way Sage carries T internally. But notice that we defined T in a way that is consistent with the text, via the use of the optional Chapter LT 100 side=’right’ keyword. If you rework examples from the text, or use Sage to assist with exercises, be sure to remember this option. In particular, when the matrix is square it might be easy to miss that you have forgotten this option. Note too that Sage uses a more general term for a linear transformation, “vector space morphism.” Just mentally translate from Sage-speak to the terms we use here in the text. If the standard way that Sage prints a linear transformation is too confusing, you can get all the relevant information with a handful of commands. T.matrix(side=’right’) [-1 0 2] [ 1 3 7] [ 1 1 1] [ 2 3 5] T.domain() Vector space of dimension 3 over Rational Field T.codomain() Vector space of dimension 4 over Rational Field So we can build a linear transformation in Sage from a matrix, as promised by Theorem MBLT. Furthermore, Theorem MLTCV says there is a matrix associated with every linear transformation. This matrix is provided in Sage by the .matrix() method, where we use the option side=’right’ to be consistent with the text. Here is Example MOLT reprised, where we define the linear transformation via a Sage symbolic function and then recover the matrix of the linear transformation. x1, x2, x3 = var(’x1, x2, x3’) outputs = [3*x1 - 2*x2 + 5*x3, x1 + x2 + x3, 9*x1 - 2*x2 + 5*x3, 4*x2 ] S_symbolic(x1, x2, x3) = outputs S = linear_transformation(QQ^3, QQ^4, S_symbolic) C = S.matrix(side=’right’); C [ 3 -2 5] [ 1 1 1] [ 9 -2 5] [ 0 4 0] x = vector(QQ, [2, -3, 3]) S(x) == C*x True LTB Linear Transformations, Bases A third way to create a linear transformation in Sage is to provide a list of images for a basis, as motivated by Theorem LTDB. The default is to use the standard basis as the inputs (Definition SUV). We will, once again, create our running example. Chapter LT 101 U = QQ^3 V = QQ^4 v1 = vector(QQ, [-1, 1, 1, 2]) v2 = vector(QQ, [ 0, 3, 1, 3]) v3 = vector(QQ, [ 2, 7, 1, 5]) T = linear_transformation(U, V, [v1, v2, v3]) T Vector space morphism represented by the matrix: [-1 1 1 2] [ 0 3 1 3] [ 2 7 1 5] Domain: Vector space of dimension 3 over Rational Field Codomain: Vector space of dimension 4 over Rational Field Notice that there is no requirement that the list of images (in Sage or in Theorem LTDB) is a basis. They do not even have to be different. They could all be the zero vector (try it). If we want to use an alternate basis for the domain, it is possible, but there are two caveats. The first caveat is that we must be sure to provide a basis for the domain, Sage will give an error if the proposed basis is not linearly independent and we are responsible for providing the right number of vectors (which should be easy). We have seen that vector spaces can have alternate bases, which prints as a “user basis.” Here will provide the domain with an alternate basis. The relevant command will create a subspace, but for now, we need to provide a big enough set to create the entire domain. It is possible to use fewer linearly independent vectors, and create a proper subspace, but then we will not be able to use this proper subspace to build the linear transformation we want. u1 = vector(QQ, [ 1, 3, -4]) u2 = vector(QQ, [-1, -2, 3]) u3 = vector(QQ, [ 1, 1, -3]) U = (QQ^3).subspace_with_basis([u1, u2, u3]) U == QQ^3 True U.basis_matrix() [ 1 3 -4] [-1 -2 3] [ 1 1 -3] U.echelonized_basis_matrix() [1 0 0] [0 1 0] [0 0 1] We can use this alternate version of U to create a linear transformation from specified images. Superficially there is nothing real special about our choices for v1, v2, v3. Chapter LT 102 V = QQ^4 v1 = vector(QQ, [-9, -18, 0, -9]) v2 = vector(QQ, [ 7, 14, 0, 7]) v3 = vector(QQ, [-7, -17, -1, -10]) Now we create the linear transformation. Here is the second caveat: the ma- trix of the linear transformation is no longer that provided by Theorem MLTCV. It may be obvious where the matrix comes from, but a full understanding of its interpretation will have to wait until Section MR. S = linear_transformation(U, V, [v1, v2, v3]) S.matrix(side=’right’) [ -9 7 -7] [-18 14 -17] [ 0 0 -1] [ -9 7 -10] We suggested our choices for v1, v2, v3 were “random.” Not so — the linear transformation S just created is equal to the linear transformation T above. If you have run all the input in this subsection, in order, then you should be able to compare the functions S and T. The next command should always produce True. u = random_vector(QQ, 3) T(u) == S(u) True Notice that T == S may not do what you expect here. Instead, the linear trans- formation method .is_equal_function() will perform a conclusive check of equal- ity of two linear transformations as functions. T.is_equal_function(S) True Can you reproduce this example? In other words, define some linear transfor- mation, any way you like. Then give the domain an alternate basis and concoct the correct images to create a second linear transformation (by the method of this subsection) which is equal to the first. PI Pre-Images Sage handles pre-images just a bit differently than our approach in the text. For the moment, we can obtain a single vector in the set that is the pre-image via the .preimage_representative() method. Understand that this method will return just one element of the pre-image set, and we have no real control over which one. Also, it is certainly possible that a pre-image is the empty set — in this case, the method will raise a ValueError. We will use our running example to illustrate. A = matrix(QQ, [[-1, 0, 2], [ 1, 3, 7], [ 1, 1, 1], [ 2, 3, 5]]) T = linear_transformation(QQ^3, QQ^4, A, side=’right’) v = vector(QQ, [1, 2, 0, 1]) u = T.preimage_representative(v) Chapter LT 103 u (-1, 1, 0) T(u) == v True T.preimage_representative(vector(QQ, [1, 2, 1, 1])) Traceback (most recent call last): ... ValueError: element is not in the image Remember, we have defined the pre-image as a set, and Sage just gives us a single element of the set. We will see in Sage ILT that the upcoming Theorem KPI explains why this is no great shortcoming in Sage. OLT Operations on Linear Transformations It is possible in Sage to add linear transformations (Definition LTA), multiply them by scalars (Definition LTSM) and compose (Definition LTC) them. Then Theorem SLTLT Theorem MLTLT, and Theorem CLTLT (respectively) tell us the results are again linear transformations. Here are some examples: U = QQ^4 V = QQ^2 A = matrix(QQ, 2, 4, [[-1, 3, 4, 5], [ 2, 0, 3, -1]]) T = linear_transformation(U, V, A, side=’right’) B = matrix(QQ, 2, 4, [[-7, 4, -2, 0], [ 1, 1, 8, -3]]) S = linear_transformation(U, V, B, side=’right’) P = S + T P Vector space morphism represented by the matrix: [-8 3] [ 7 1] [ 2 11] [ 5 -4] Domain: Vector space of dimension 4 over Rational Field Codomain: Vector space of dimension 2 over Rational Field Q = S*5 Q Vector space morphism represented by the matrix: [-35 5] [ 20 5] [-10 40] [ 0 -15] Domain: Vector space of dimension 4 over Rational Field Chapter LT 104 Codomain: Vector space of dimension 2 over Rational Field Perhaps the only surprise in all this is the necessity of writing scalar multipli- cation on the right of the linear transformation (rather on the left, as we do in the text). We will recycle the linear transformation T from above and redefine S to form an example of composition. W = QQ^3 C = matrix(QQ, [[ 4, -2], [-1, 3], [-3, 2]]) S = linear_transformation(V, W, C, side=’right’) R = S*T R Vector space morphism represented by the matrix: [ -8 7 7] [ 12 -3 -9] [ 10 5 -6] [ 22 -8 -17] Domain: Vector space of dimension 4 over Rational Field Codomain: Vector space of dimension 3 over Rational Field We use the star symbol (*) to indicate composition of linear transformations. Notice that the order of the two linear transformations we compose is important, and Sage’s order agrees with the text. The order does not have to agree, and there are good arguments to have it reversed, so be careful if you read about composition elsewhere. This is a good place to expand on Theorem VSLT, which says that with defini- tions of addition and scalar multiplication of linear transformations we then arrive at a vector space. A vector space full of linear transformations. Objects in Sage have “parents” — vectors have vector spaces for parents, fractions of integers have the rationals as parents. What is the parent of a linear transformation? Let us see, by investigating the parent of S just defined above. P = S.parent() P Set of Morphisms (Linear Transformations) from Vector space of dimension 2 over Rational Field to Vector space of dimension 3 over Rational Field “Morphism” is a general term for a function that “preserves structure” or “re- spects operations.” In Sage a collection of morphisms is referenced as a “homset” or a “homspace.” In this example, we have a homset that is the vector space of linear transformations that go from a dimension 2 vector space over the rationals to a dimension 3 vector space over the rationals. What can we do with it? A few things, but not everything you might imagine. It does have a basis, containing a few very simple linear transformations: P.basis() (Vector space morphism represented by the matrix: [1 0 0] [0 0 0] Domain: Vector space of dimension 2 over Rational Field Codomain: Vector space of dimension 3 over Rational Field, Vector space morphism represented by the matrix: Chapter LT 105 [0 1 0] [0 0 0] Domain: Vector space of dimension 2 over Rational Field Codomain: Vector space of dimension 3 over Rational Field, Vector space morphism represented by the matrix: [0 0 1] [0 0 0] Domain: Vector space of dimension 2 over Rational Field Codomain: Vector space of dimension 3 over Rational Field, Vector space morphism represented by the matrix: [0 0 0] [1 0 0] Domain: Vector space of dimension 2 over Rational Field Codomain: Vector space of dimension 3 over Rational Field, Vector space morphism represented by the matrix: [0 0 0] [0 1 0] Domain: Vector space of dimension 2 over Rational Field Codomain: Vector space of dimension 3 over Rational Field, Vector space morphism represented by the matrix: [0 0 0] [0 0 1] Domain: Vector space of dimension 2 over Rational Field Codomain: Vector space of dimension 3 over Rational Field) You can create a set of linear transformations with the Hom() function, simply by giving the domain and codomain. H = Hom(QQ^6, QQ^9) H Set of Morphisms (Linear Transformations) from Vector space of dimension 6 over Rational Field to Vector space of dimension 9 over Rational Field An understanding of Sage’s homsets is not critical to understanding the use of Sage during the remainder of this course. But such an understanding can be very useful in understanding some of Sage’s more advanced and powerful features. Section ILT: Injective Linear Transformations ILT Injective Linear Transformations By now, you have probably already figured out how to determine if a linear trans- formation is injective, and what its kernel is. You may also now begin to understand why Sage calls the null space of a matrix a kernel. Here are two examples, first a reprise of Example NKAO. U = QQ^3 V = QQ^5 x1, x2, x3 = var(’x1, x2, x3’) outputs = [ -x1 + x2 - 3*x3, -x1 + 2*x2 - 4*x3, x1 + x2 + x3, 2*x1 + 3*x2 + x3, x1 + 2*x3] T_symbolic(x1, x2, x3) = outputs T = linear_transformation(U, V, T_symbolic) Chapter LT 106 T.is_injective() False T.kernel() Vector space of degree 3 and dimension 1 over Rational Field Basis matrix: [ 1 -1/2 -1/2] So we have a concrete demonstration of one half of Theorem KILT. Here is the second example, a do-over for Example TKAP, but renamed as S. U = QQ^3 V = QQ^5 x1, x2, x3 = var(’x1, x2, x3’) outputs = [ -x1 + x2 + x3, -x1 + 2*x2 + 2*x3, x1 + x2 + 3*x3, 2*x1 + 3*x2 + x3, -2*x1 + x2 + 3*x3] S_symbolic(x1, x2, x3) = outputs S = linear_transformation(U, V, S_symbolic) S.is_injective() True S.kernel() Vector space of degree 3 and dimension 0 over Rational Field Basis matrix: [] S.kernel() == U.subspace([]) True And so we have a concrete demonstration of the other half of Theorem KILT. Now that we have Theorem KPI, we can return to our discussion from Sage PI. The .preimage_representative() method of a linear transformation will give us a single element of the pre-image, with no other guarantee about the nature of that element. That is fine, since this is all Theorem KPI requires (in addition to the kernel). Remember that not every element of the codomain may have a nonempty pre-image (as indicated in the hypotheses of Theorem KPI). Here is an example using T from above, with a choice of a codomain element that has a nonempty pre-image. TK = T.kernel() v = vector(QQ, [2, 3, 0, 1, -1]) u = T.preimage_representative(v) u Chapter LT 107 (-1, 1, 0) Now the following will create random elements of the preimage of v, which can be verified by the test always returning True. Use the compute cell just below if you are curious what p looks like. p = u + TK.random_element() T(p) == v True p # random (-13/10, 23/20, 3/20) As suggested, some choices of v can lead to empty pre-images, in which case Theorem KPI does not even apply. v = vector(QQ, [4, 6, 3, 1, -2]) u = T.preimage_representative(v) Traceback (most recent call last): ... ValueError: element is not in the image The situation is less interesting for an injective linear transformation. Still, pre- images may be empty, but when they are nonempty, they are just singletons (a single element) since the kernel is empty. So a repeat of the above example, with S rather than T, would not be very informative. CILT Composition of Injective Linear Transformations One way to use Sage is to construct examples of theorems and verify the conclusions. Sometimes you will get this wrong: you might build an example that does not satisfy the hypotheses, or your example may not satisfy the conclusions. This may be because you are not using Sage properly, or because you do not understand a definition or a theorem, or in very limited cases you may have uncovered a bug in Sage (which is always the preferred explanation!). But in the process of trying to understand a discrepancy or unexpected result, you will learn much more, both about linear algebra and about Sage. And Sage is incredibly patient — it will stay up with you all night to help you through a rough patch. Let us illustrate the above in the context of Theorem CILTI. The hypotheses indicate we need two injective linear transformations. Where will get two such linear transformations? Well, the contrapositive of Theorem ILTD tells us that if the dimension of the domain exceeds the dimension of the codomain, we will never be injective. So we should at a minimum avoid this scenario. We can build two linear transformations from matrices created randomly, and just hope that they lead to injective linear transformations. Here is an example of how we create examples like this. The random matrix has single-digit entries, and almost always will lead to an injective linear transformation, though we cannot be absolutely certain. Evaluate this cell repeatedly, to see how rarely the result is not injective. E = random_matrix(ZZ, 3, 2, x=-9, y=9) T = linear_transformation(QQ^2, QQ^3, E, side=’right’) T.is_injective() # random Chapter LT 108 True Our concrete example below was created this way, so here we go. U = QQ^2 V = QQ^3 W = QQ^4 A = matrix(QQ, 3, 2, [[-4, -1], [-3, 3], [ 7, -6]]) B = matrix(QQ, 4, 3, [[ 7, 0, -1], [ 6, 2, 7], [ 3, -1, 2], [-6, -1, 1]]) T = linear_transformation(U, V, A, side=’right’) T.is_injective() True S = linear_transformation(V, W, B, side=’right’) S.is_injective() True C = S*T C.is_injective() True Section SLT: Surjective Linear Transformations SLT Surjective Linear Transformations The situation in Sage for surjective linear transformations is similar to that for injective linear transformations. One distinction — what your text calls the range of a linear transformation is called the image of a transformation, obtained with the .image() method. Sage’s term is more commonly used, so you are likely to see it again. As before, two examples, first up is Example RAO. U = QQ^3 V = QQ^5 x1, x2, x3 = var(’x1, x2, x3’) outputs = [ -x1 + x2 - 3*x3, -x1 + 2*x2 - 4*x3, x1 + x2 + x3, 2*x1 + 3*x2 + x3, x1 + 2*x3] T_symbolic(x1, x2, x3) = outputs T = linear_transformation(U, V, T_symbolic) T.is_surjective() False Chapter LT 109 T.image() Vector space of degree 5 and dimension 2 over Rational Field Basis matrix: [ 1 0 -3 -7 -2] [ 0 1 2 5 1] Besides showing the relevant commands in action, this demonstrates one half of Theorem RSLT. Now a reprise of Example FRAN. U = QQ^5 V = QQ^3 x1, x2, x3, x4, x5 = var(’x1, x2, x3, x4, x5’) outputs = [2*x1 + x2 + 3*x3 - 4*x4 + 5*x5, x1 - 2*x2 + 3*x3 - 9*x4 + 3*x5, 3*x1 + 4*x3 - 6*x4 + 5*x5] S_symbolic(x1, x2, x3, x4, x5) = outputs S = linear_transformation(U, V, S_symbolic) S.is_surjective() True S.image() Vector space of degree 3 and dimension 3 over Rational Field Basis matrix: [1 0 0] [0 1 0] [0 0 1] S.image() == V True Previously, we have chosen elements of the codomain which have nonempty or empty preimages. We can now explain how to do this predictably. Theorem RPI explains that elements of the codomain with nonempty pre-images are precisely elements of the image. Consider the non-surjective linear transformation T from above. TI = T.image() B = TI.basis() B [ (1, 0, -3, -7, -2), (0, 1, 2, 5, 1) ] b0 = B[0] b1 = B[1] Chapter LT 110 Now any linear combination of the basis vectors b0 and b1 must be an element of the image. Moreover, the first two slots of the resulting vector will equal the two scalars we choose to create the linear combination. But most importantly, see that the three remaining slots will be uniquely determined by these two choices. This means there is exactly one vector in the image with these values in the first two slots. So if we construct a new vector with these two values in the first two slots, and make any part of the last three slots even slightly different, we will form a vector that cannot be in the image, and will thus have a preimage that is an empty set. Whew, that is probably worth reading carefully several times, perhaps in conjunction with the example following. in_image = 4*b0 + (-5)*b1 in_image (4, -5, -22, -53, -13) T.preimage_representative(in_image) (-13, -9, 0) outside_image = 4*b0 + (-5)*b1 + vector(QQ, [0, 0, 0, 0, 1]) outside_image (4, -5, -22, -53, -12) T.preimage_representative(outside_image) Traceback (most recent call last): ... ValueError: element is not in the image CSLT Composition of Surjective Linear Transformations As we mentioned in the last section, experimenting with Sage is a worthwhile com- plement to other methods of learning mathematics. We have purposely avoided pro- viding illustrations of deeper results, such as Theorem ILTB and Theorem SLTB, which you should now be equipped to investigate yourself. For completeness, and since composition will be very important in the next few sections, we will provide an illustration of Theorem CSLTS. Similar to what we did in the previous sec- tion, we choose dimensions suggested by Theorem SLTD, and then use randomly constructed matrices to form a pair of surjective linear transformations. U = QQ^4 V = QQ^3 W = QQ^2 A = matrix(QQ, 3, 4, [[ 3, -2, 8, -9], [-1, 3, -4, -1], [ 3, 2, 8, 3]]) T = linear_transformation(U, V, A, side=’right’) T.is_surjective() Chapter LT 111 True B = matrix(QQ, 2, 3, [[ 8, -5, 3], [-2, 1, 1]]) S = linear_transformation(V, W, B, side=’right’) S.is_surjective() True C = S*T C.is_surjective() True Section IVLT: Invertible Linear Transformations IVLT Invertible Linear Transformations Of course, Sage can compute the inverse of a linear transformation. However, not every linear transformation has an inverse, and we will see shortly how to determine this. For now, take this example as just an illustration of the basics (both mathematically and for Sage). U = QQ^4 V = QQ^4 x1, x2, x3, x4 = var(’x1, x2, x3, x4’) outputs = [ x1 + 2*x2 - 5*x3 - 7*x4, x2 - 3*x3 - 5*x4, x1 + 2*x2 - 4*x3 - 6*x4, -2*x1 - 2*x2 + 7*x3 + 8*x4 ] T_symbolic(x1, x2, x3, x4) = outputs T = linear_transformation(U, V, T_symbolic) S = T.inverse() S Vector space morphism represented by the matrix: [-8 7 -6 5] [ 2 -3 2 -2] [ 5 -3 4 -3] [-2 2 -1 1] Domain: Vector space of dimension 4 over Rational Field Codomain: Vector space of dimension 4 over Rational Field We can build the composition of T and its inverse, S, in both orders. We will optimistically name these as identity linear transformations, as predicted by Defi- nition IVLT. Run the cells to define the compositions, then run the compute cells with the random vectors repeatedly — they should always return True. IU = S*T IV = T*S u = random_vector(QQ, 4) IU(u) == u # random Chapter LT 112 True v = random_vector(QQ, 4) IV(v) == v # random True We can also check that the compositions are the same as the identity linear transformation itself. We will do one, you can try the other. id = linear_transformation(U, U, identity_matrix(QQ, 4)) IU.is_equal_function(id) True CIVLT Computing the Inverse of a Linear Transformations Theorem ILTIS gives us a straightforward condition equivalence for an invertible linear transformation, but of course, it is even easier in Sage. U = QQ^4 V = QQ^4 x1, x2, x3, x4 = var(’x1, x2, x3, x4’) outputs = [ x1 + 2*x2 - 5*x3 - 7*x4, x2 - 3*x3 - 5*x4, x1 + 2*x2 - 4*x3 - 6*x4, -2*x1 - 2*x2 + 7*x3 + 8*x4 ] T_symbolic(x1, x2, x3, x4) = outputs T = linear_transformation(U, V, T_symbolic) T.is_invertible() True As easy as that is, it is still instructive to walk through an example similar to Example CIVLT using Sage, as a further illustration of the second half of the proof of Theorem ILTIS. Since T is surjective, every element of the codomain has a nonempty pre-image, and since T is injective, the pre-image of each element is a single element. Keep these facts in mind and convince yourself that the procedure below would never raise an error, and always has a unique result. We first compute the pre-image of each element of a basis of the codomain. preimages = [T.preimage_representative(v) for v in V.basis()] preimages [(-8, 7, -6, 5), (2, -3, 2, -2), (5, -3, 4, -3), (-2, 2, -1, 1)] Then we define a new linear transformation, from V to U, which turns it around and uses the preimages as a set of images defining the new linear transformation. Explain to yourself how we know that preimages is a basis for U, and why this will create an invertible linear transformation. S = linear_transformation(V, U, preimages) S Chapter LT 113 Vector space morphism represented by the matrix: [-8 7 -6 5] [ 2 -3 2 -2] [ 5 -3 4 -3] [-2 2 -1 1] Domain: Vector space of dimension 4 over Rational Field Codomain: Vector space of dimension 4 over Rational Field S.is_equal_function(T.inverse()) True While this is a simple two-step procedure (form preimages, construct linear transformation), realize that this is not the process that Sage uses internally. Notice that the essence of this construction is that when we work with an invert- ible linear transformation, the method .preimage_representative() behaves as a function (we mean the precise mathematical definition here) — it is always defined and always produces just one well-defined output. Here the linear_transformation() constructor is extending it to a linear function based on its action on a (finite) basis of the domain. LTOE Linear Transformation Odds and Ends We should mention that the notation T^-1 will yield an inverse of a linear transfor- mation in Sage. U = QQ^4 V = QQ^4 x1, x2, x3, x4 = var(’x1, x2, x3, x4’) outputs = [ x1 + 2*x2 - 5*x3 - 7*x4, x2 - 3*x3 - 5*x4, x1 + 2*x2 - 4*x3 - 6*x4, -2*x1 - 2*x2 + 7*x3 + 8*x4] T_symbolic(x1, x2, x3, x4) = outputs T = linear_transformation(U, V, T_symbolic) T^-1 Vector space morphism represented by the matrix: [-8 7 -6 5] [ 2 -3 2 -2] [ 5 -3 4 -3] [-2 2 -1 1] Domain: Vector space of dimension 4 over Rational Field Codomain: Vector space of dimension 4 over Rational Field Also, the rank and nullity are what you might expect. Recall that for a matrix Sage provides a left nullity and a right nullity. There is no such distinction for linear transformations. We verify Theorem RPNDD as an example. U = QQ^3 V = QQ^4 x1, x2, x3 = var(’x1, x2, x3’) outputs = [ -x1 + 2*x3, x1 + 3*x2 + 7*x3, x1 + x2 + x3, 2*x1 + 3*x2 + 5*x3] R_symbolic(x1, x2, x3) = outputs Chapter LT 114 R = linear_transformation(QQ^3, QQ^4, R_symbolic) R.rank() 2 R.nullity() 1 R.rank() + R.nullity() == U.dimension() True SUTH1 Sage Under The Hood, Round 1 We can parallel the above discussion about systems of linear equations and linear transformations using Sage. We begin with a matrix that we will use as a coef- ficient matrix for systems of equations, and then use the same matrix to define the associated linear transformation (acting on vectors placed to the right of the matrix). A = matrix(QQ, [[-1, 0, 2], [ 1, 3, 7], [ 1, 1, 1], [ 2, 3, 5]]) T = linear_transformation(QQ^3, QQ^4, A, side=’right’) We solve a linear system using the coefficient matrix, and compute an element of the pre-image of the linear transformation. v = vector(QQ, [1, 2, 0, 1]) A.solve_right(v) (-1, 1, 0) T.preimage_representative(v) (-1, 1, 0) We compute a null space of the coefficient matrix and a kernel of the linear transformation, so as to understand the full solution set or the full preimage set. A.right_kernel() Vector space of degree 3 and dimension 1 over Rational Field Basis matrix: [ 1 -3/2 1/2] T.kernel() Chapter LT 115 Vector space of degree 3 and dimension 1 over Rational Field Basis matrix: [ 1 -3/2 1/2] We compute a column space of the coefficient matrix and an image (range) of the linear transformation to help us build an inconsistent system. A.column_space() Vector space of degree 4 and dimension 2 over Rational Field Basis matrix: [ 1 0 -2/3 -1] [ 0 1 1/3 1] T.image() Vector space of degree 4 and dimension 2 over Rational Field Basis matrix: [ 1 0 -2/3 -1] [ 0 1 1/3 1] A vector that creates a system with no solution is a vector that has an empty pre-image. v = vector(QQ, [1, 2, 1, 1]) A.solve_right(v) Traceback (most recent call last): ... ValueError: matrix equation has no solutions T.preimage_representative(v) Traceback (most recent call last): ... ValueError: element is not in the image Note that we could redo the above, replacing uses of “right” by “left” and uses of “column” by “row.” The output would not change. We suggest in the text that one could develop the theory of linear transforma- tions from scratch, and then obtain many of our initial results about systems of equations and matrices as a consequence. In Sage it is the reverse. Sage imple- ments many advanced concepts from various areas of mathematics by translating fundamental computations into the language of linear algebra. In turn, many of Sage’s linear algebra routines ultimately depend on very fast algorithms for basic operations on matrices, such as computing an echelon form, a null space, or a span. Chapter R Representations Section VR: Vector Representations VR Vector Representations Vector representation is described in the text in a fairly abstract fashion. Sage will support this view (which will be useful in the next section), as well as provide a more practical approach. We will explain both approaches. We begin with an arbitrarily chosen basis. We then create an alternate version of QQ^4 with this basis as a“user basis”, namely V. v0 = vector(QQ, [ 1, 1, 1, 0]) v1 = vector(QQ, [ 1, 2, 3, 2]) v2 = vector(QQ, [ 2, 2, 3, 2]) v3 = vector(QQ, [-1, 3, 5, 5]) B = [v0, v1, v2, v3] V = (QQ^4).subspace_with_basis(B) V Vector space of degree 4 and dimension 4 over Rational Field User basis matrix: [ 1 1 1 0] [ 1 2 3 2] [ 2 2 3 2] [-1 3 5 5] V.echelonized_basis_matrix() [1 0 0 0] [0 1 0 0] [0 0 1 0] [0 0 0 1] Now, the construction of a linear transformation will use the basis provided for V. In the proof of Theorem VRLT we defined a linear transformation T that equaled ρB . T was defined by taking the basis vectors of B to the basis composed of standard unit vectors (Definition SUV). This is exactly what we will accomplish in the following construction. Note how the basis associated with the domain is automatically paired with the elements of the basis for the codomain. rho = linear_transformation(V, QQ^4, (QQ^4).basis()) rho 116 Chapter R 117 Vector space morphism represented by the matrix: [1 0 0 0] [0 1 0 0] [0 0 1 0] [0 0 0 1] Domain: Vector space of degree 4 and dimension 4 over Rational Field User basis matrix: [ 1 1 1 0] [ 1 2 3 2] [ 2 2 3 2] [-1 3 5 5] Codomain: Vector space of dimension 4 over Rational Field First, we verify Theorem VRILT: rho.is_invertible() True Notice that the matrix of the linear transformation is the identity matrix. This might look odd now, but we will have a full explanation soon. Let us see if this linear transformation behaves as it should. We will “coordinatize” an arbitrary vector, w. w = vector(QQ, [-13, 28, 45, 43]) rho(w) (3, 5, -6, 9) lincombo = 3*v0 + 5*v1 + (-6)*v2 + 9*v3 lincombo (-13, 28, 45, 43) lincombo == w True Notice how the expression for lincombo is exactly the messy expression displayed in Definition VR. More precisely, we could even write this as: w == sum([rho(w)[i]*B[i] for i in range(4)]) True Or we can test this equality repeatedly with random vectors. u = random_vector(QQ, 4) u == sum([rho(u)[i]*B[i] for i in range(4)]) True Finding a vector representation is such a fundamental operation that Sage has an easier command, bypassing the need to create a linear transformation. It does still require constructing a vector space with the alternate basis. Here goes, repeating the prior example. Chapter R 118 w = vector(QQ, [-13, 28, 45, 43]) V.coordinate_vector(w) (3, 5, -6, 9) Boom! SUTH2 Sage Under The Hood, Round 2 You will have noticed that we have never constructed examples involving our favorite abstract vector spaces, such as the vector space of polynomials with fixed maximum degree, the vector space of matrices of a fixed size, or even the crazy vector space. There is nothing to stop us (or you) from implementing these examples in Sage as vector spaces. Maybe someday it will happen. But since Sage is built to be a tool for serious mathematical research, the designers recognize that this is not necessary. Theorem CFDVS tells us that every finite-dimensional vector space can be de- scribed (loosely speaking) by just a field of scalars (for us, C in the text, QQ in Sage) and the dimension. You can study whatever whacky vector space you might dream up, or whatever very complicated vector space that is important for particle physics, and through vector representation (“coordinatization”), you can convert your calculations to and from Sage. Section MR: Matrix Representations MR Matrix Representations It is very easy to create matrix representations of linear transformations in Sage. Our examples in the text have used abstract vector spaces to make it a little easier to keep track of the domain and codomain. With Sage we will have to consistently use variants of QQ^n, but we will use non-obvious bases to make it nontrivial and interesting. Here are the components of our first example, one linear function, two vector spaces, four bases. x1, x2, x3, x4 = var(’x1, x2, x3, x4’) outputs = [3*x1 + 7*x2 + x3 - 4*x4, 2*x1 + 5*x2 + x3 - 3*x4, -x1 - 2*x2 + x4] T_symbolic(x1, x2, x3, x4) = outputs U = QQ^4 V = QQ^3 b0 = vector(QQ, [ 1, 1, -1, 0]) b1 = vector(QQ, [-1, 0, -2, 7]) b2 = vector(QQ, [ 0, 1, -2, 4]) b3 = vector(QQ, [-2, 0, -1, 6]) B = [b0, b1, b2, b3] c0 = vector(QQ, [ 1, 6, -6]) c1 = vector(QQ, [ 0, 1, -1]) c2 = vector(QQ, [-2, -3, 4]) C = [c0, c1, c2] d0 = vector(QQ, [ 1, -3, 2, -1]) d1 = vector(QQ, [ 0, 1, 0, 1]) d2 = vector(QQ, [-1, 2, -1, -1]) d3 = vector(QQ, [ 2, -8, 4, -3]) D = [d0, d1, d2, d3] e0 = vector(QQ, [ 0, 1, -3]) e1 = vector(QQ, [-1, 2, -1]) e2 = vector(QQ, [ 2, -4, 3]) Chapter R 119 E = [e0, e1, e2] We create alternate versions of the domain and codomain by providing them with bases other than the standard ones. Then we build the linear transformation and ask for its .matrix(). We will use numerals to distinguish our two examples. U1 = U.subspace_with_basis(B) V1 = V.subspace_with_basis(C) T1 = linear_transformation(U1, V1, T_symbolic) T1.matrix(side=’right’) [ 15 -67 -25 -61] [-75 326 120 298] [ 3 -17 -7 -15] Now, we do it again, but with the bases D and E. U2 = U.subspace_with_basis(D) V2 = V.subspace_with_basis(E) T2 = linear_transformation(U2, V2, T_symbolic) T2.matrix(side=’right’) [ -32 8 38 -91] [-148 37 178 -422] [ -80 20 96 -228] So once the pieces are in place, it is quite easy to obtain a matrix representation. You might experiment with examples from the text, by first mentally converting elements of abstract vector spaces into column vectors via vector representations using simple bases for the abstract spaces. The linear transformations T1 and T2 are not different as functions, despite the fact that Sage treats them as unequal (since they have different matrix representa- tions). We can check that the two functions behave identically, first with random testing. Repeated execution of the following compute cell should always produce True. u = random_vector(QQ, 4) T1(u) == T2(u) True A better approach would be to see if T1 and T2 agree on a basis for the domain, and to avoid any favoritism, we will use the basis of U itself. Convince yourself that this is a proper application of Theorem LTDB to demonstrate equality. all([T1(u) == T2(u) for u in U.basis()]) True Or we can just ask if they are equal functions (a method that is implemented using the previous idea about agreeing on a basis). T1.is_equal_function(T2) True We can also demonstrate Theorem FTMR — we will use the representation for T1. We need a couple of vector representation linear transformations. Chapter R 120 rhoB = linear_transformation(U1, U, U.basis()) rhoC = linear_transformation(V1, V, V.basis()) Now we experiment with a “random” element of the domain. Study the third computation carefully, as it is the Sage version of Theorem FTMR. You might compute some of the pieces of this expression to build up the final expression, and you might duplicate the experiment with a different input and/or with the T2 version. T_symbolic(-3, 4, -5, 2) (6, 3, -3) u = vector(QQ, [-3, 4, -5, 2]) T1(u) (6, 3, -3) (rhoC^-1)(T1.matrix(side=’right’)*rhoB(u)) (6, 3, -3) One more time, but we will replace the vector representation linear transforma- tions with Sage conveniences. Recognize that the .linear_combination_of_basis() method is the inverse of vector representation (it is “un-coordinatization”). output_coord = T1.matrix(side=’right’)*U1.coordinate_vector(u) V1.linear_combination_of_basis(output_coord) (6, 3, -3) We have concentrated here on the first version of the conclusion of Theorem FTMR. You could experiment with the second version in a similar manner. Extra credit: what changes do you need to make in any of these experiments if you remove the side=’right’ option? SUTH3 Sage Under The Hood, Round 3 We have seen that Sage is able to add two linear transformations, or multiply a single linear transformation by a scalar. Under the hood, this is accomplished by simply adding matrix representations, or multiplying a matrix by a scalar, according to Theorem MRSLT and Theorem MRMLT (respectively). Theorem MRCLT says linear transformation composition is matrix representation multiplication. Is it still a mystery why we use the symbol * for linear transformation composition in Sage? We could do several examples here, but you should now be able to construct them yourselves. We will do just one, linear transformation composition is matrix representation multiplication. x1, x2, x3, x4 = var(’x1, x2, x3, x4’) outputs = [-5*x1 - 2*x2 + x3, 4*x1 - 3*x2 - 3*x3, 4*x1 - 6*x2 - 4*x3, 5*x1 + 3*x2 ] T_symbolic(x1, x2, x3) = outputs outputs = [-3*x1 - x2 + x3 + 2*x4, 7*x1 + x2 + x3 - x4] Chapter R 121 S_symbolic(x1, x2, x3, x4) = outputs b0 = vector(QQ, [-1, -2, 2]) b1 = vector(QQ, [ 1, 1, 0]) b2 = vector(QQ, [ 0, 3, -5]) U = (QQ^3).subspace_with_basis([b0, b1, b2]) c0 = vector(QQ, [ 0, 0, 2, 1]) c1 = vector(QQ, [ 2, -3, -1, -6]) c2 = vector(QQ, [-2, 3, 2, 7]) c3 = vector(QQ, [ 1, -2, -4, -6]) V = (QQ^4).subspace_with_basis([c0, c1, c2, c3]) d0 = vector(QQ, [3, 4]) d1 = vector(QQ, [2, 3]) W = (QQ^2).subspace_with_basis([d0, d1]) T = linear_transformation(U, V, T_symbolic) S = linear_transformation(V, W, S_symbolic) (S*T).matrix(’right’) [-321 218 297] [ 456 -310 -422] S.matrix(side=’right’)*T.matrix(side=’right’) [-321 218 297] [ 456 -310 -422] Extra credit: what changes do you need to make if you dropped the side=’right’ option on these three matrix representations? LTR Linear Transformation Restrictions Theorem KNSI and Theorem RCSI have two of the most subtle proofs we have seen so far. The conclusion that two vector spaces are isomorphic is established by actually constructing an isomorphism between the vector spaces. To build the isomorphism, we begin with a familar object, a vector representation linear trans- formation, but the hard work is showing that we can “restrict” the domain and codomain of this function and still arrive at a legitimate (invertible) linear trans- formation. In an effort to make the proofs more concrete, we will walk through a nontrivial example for Theorem KNSI, and you might try to do the same for The- orem RCSI. (An understanding of this subsection is not needed for the remainder — its second purpose is to demonstrate some of the powerful tools Sage provides.) Here are the pieces. We build a linear transformation with two different repre- sentations, one with respect to standard bases, the other with respect to less-obvious bases. x1, x2, x3, x4, x5 = var(’x1, x2, x3, x4, x5’) outputs = [ x1 - x2 - 5*x3 + x4 + x5, x1 - 2*x3 - x4 - x5, - x2 - 3*x3 + 2*x4 + 2*x5, -x1 + x2 + 5*x3 - x4 - x5] T_symbolic(x1, x2, x3, x4, x5) = outputs b0 = vector(QQ, [-1, 6, 5, 5, 1]) b1 = vector(QQ, [-1, 5, 4, 4, 1]) b2 = vector(QQ, [-2, 4, 3, 2, 5]) b3 = vector(QQ, [ 1, -1, 0, 1, -5]) b4 = vector(QQ, [ 3, -7, -6, -5, -4]) U = (QQ^5).subspace_with_basis([b0, b1, b2, b3, b4]) c0 = vector(QQ, [1, 1, 1, -3]) Chapter R 122 c1 = vector(QQ, [-2, 3, -6, -7]) c2 = vector(QQ, [0, -1, 1, 2]) c3 = vector(QQ, [-1, 3, -4, -7]) V = (QQ^4).subspace_with_basis([c0, c1, c2, c3]) T_plain = linear_transformation(QQ^5, QQ^4, T_symbolic) T_fancy = linear_transformation( U, V, T_symbolic) Now we compute the kernel of the linear transformation using the “plain” ver- sion, and the null space of a matrix representation coming from the “fancy” version. K = T_plain.kernel() K Vector space of degree 5 and dimension 3 over Rational Field Basis matrix: [ 1 0 2/7 0 3/7] [ 0 1 -1/7 0 2/7] [ 0 0 0 1 -1] MK = T_fancy.matrix(side=’right’).right_kernel() MK Vector space of degree 5 and dimension 3 over Rational Field Basis matrix: [ 1 0 0 -97/203 164/203] [ 0 1 0 -10/29 19/29] [ 0 0 1 129/203 100/203] So quite obviously, the kernel of the linear transformation is quite different look- ing from the null space of the matrix representation. Though it is no accident that they have the same dimension. Now we build the necessary vector representation, and use two Sage commands to “restrict” the function to a smaller domain (the kernel of the linear transformation) and a smaller codomain (the null space of the matrix representation relative to nonstandard bases). rhoB = linear_transformation(U, QQ^5, (QQ^5).basis()) rho_restrict = rhoB.restrict_domain(K).restrict_codomain(MK) rho_restrict Vector space morphism represented by the matrix: [ 33/7 -37/7 -11/7] [-13/7 22/7 -12/7] [ -4 5 6] Domain: Vector space of degree 5 and dimension 3 over Rational Field Basis matrix: [ 1 0 2/7 0 3/7] [ 0 1 -1/7 0 2/7] [ 0 0 0 1 -1] Codomain: Vector space of degree 5 and dimension 3 over Rational Field Basis matrix: [ 1 0 0 -97/203 164/203] [ 0 1 0 -10/29 19/29] [ 0 0 1 129/203 100/203] The first success is that the restriction was even created. Sage would recognize if the original linear transformation ever carried an input from the restricted domain to an output that was not contained in the proposed codomain, and would have raised Chapter R 123 an error in that event. Phew! Guaranteeing this success was the first big step in the proof of Theorem KNSI. Notice that the matrix representation of the restriction is a 3 × 3 matrix, since the restriction runs between a domain and codomain that each have dimension 3. These two vector spaces (the domain and codomain of the restriction) have dimension 3 but still contain vectors with 5 entries in their un-coordinatized versions. The next two steps of the proof show that the restriction is injective (easy in the proof) and surjective (hard in the proof). In Sage, here is the second success, rho_restrict.is_injective() True rho_restrict.is_surjective() True Verified as invertible, rho_restrict qualifies as an isomorphism between the linear transformation kernel, K, and the matrix representation null space, MK. Only an example, but still very nice. Your turn — can you create a verfication of The- orem RCSI (for this example, or some other nontrivial example you might create yourself)? NME9 Nonsingular Matrix Equivalences, Round 9 Our final fact about nonsingular matrices expresses the correspondence between invertible matrices and invertible linear transformations. As a Sage demonstration, we will begin with an invertible linear transformation and examine two matrix representations. We will create the linear transformation with nonstandard bases and then compute its representation relative to standard bases. x1, x2, x3, x4 = var(’x1, x2, x3, x4’) outputs = [ x1 - 2*x3 - 4*x4, x2 - x3 - 5*x4, -x1 - 2*x2 + 2*x3 + 7*x4, - x2 + x4] T_symbolic(x1, x2, x3, x4) = outputs b0 = vector(QQ, [ 1, -2, -1, 8]) b1 = vector(QQ, [ 0, 1, 0, -2]) b2 = vector(QQ, [-1, -2, 2, -5]) b3 = vector(QQ, [-1, -3, 2, -2]) U = (QQ^4).subspace_with_basis([b0, b1, b2, b3]) c0 = vector(QQ, [ 3, -1, 4, -8]) c1 = vector(QQ, [ 1, 0, 1, -1]) c2 = vector(QQ, [ 0, 2, -1, 6]) c3 = vector(QQ, [-1, 2, -2, 8]) V = (QQ^4).subspace_with_basis([c0, c1, c2, c3]) T = linear_transformation(U, V, T_symbolic) T Vector space morphism represented by the matrix: [ 131 -56 -321 366] [ -37 17 89 -102] [ -61 25 153 -173] [ -7 -1 24 -25] Domain: Vector space of degree 4 and dimension 4 over Rational Field User basis matrix: Chapter R 124 [ 1 -2 -1 8] [ 0 1 0 -2] [-1 -2 2 -5] [-1 -3 2 -2] Codomain: Vector space of degree 4 and dimension 4 over Rational Field User basis matrix: [ 3 -1 4 -8] [ 1 0 1 -1] [ 0 2 -1 6] [-1 2 -2 8] T.is_invertible() True T.matrix(side=’right’).is_invertible() True (T^-1).matrix(side=’right’) == (T.matrix(side=’right’))^-1 True We can convert T to a new representation using standard bases for QQ^4 by computing images of the standard basis. images = [T(u) for u in (QQ^4).basis()] T_standard = linear_transformation(QQ^4, QQ^4, images) T_standard Vector space morphism represented by the matrix: [ 1 0 -1 0] [ 0 1 -2 -1] [-2 -1 2 0] [-4 -5 7 1] Domain: Vector space of dimension 4 over Rational Field Codomain: Vector space of dimension 4 over Rational Field T_standard.matrix(side=’right’).is_invertible() True Understand that any matrix representation of T_symbolic will have an invert- ible matrix representation, no matter which bases are used. If you look at the matrix representation of T_standard and the definition of T_symbolic the con- struction of this example will be transparent, especially if you know the random matrix constructor, A = random_matrix(QQ, 4, algorithm=’unimodular’, upper_bound=9) A # random Chapter R 125 [-1 -1 2 1] [ 1 1 -1 0] [-2 -1 5 6] [ 1 0 -4 -5] Section CB: Change of Basis ENDO Endomorphisms An endomorphism is an “operation-preserving” function (a “morphism”) whose domain and codomain are equal. Sage takes this definition one step further for linear transformations and requires that the domain and codomain have the same bases (either a default echelonized basis or the same user basis). When a linear transfor- mation meets this extra requirement, several natural methods become available. Principally, we can compute the eigenvalues provided by Definition EELT. We also get a natural notion of a characteristic polynomial. x1, x2, x3, x4 = var(’x1, x2, x3, x4’) outputs = [ 4*x1 + 2*x2 - x3 + 8*x4, 3*x1 - 5*x2 - 9*x3 , 6*x2 + 7*x3 + 6*x4, -3*x1 + 2*x2 + 5*x3 - 3*x4] T_symbolic(x1, x2, x3, x4) = outputs T = linear_transformation(QQ^4, QQ^4, T_symbolic) T.eigenvalues() [3, -2, 1, 1] cp = T.characteristic_polynomial() cp x^4 - 3*x^3 - 3*x^2 + 11*x - 6 cp.factor() (x - 3) * (x + 2) * (x - 1)^2 Now the question of eigenvalues being elements of the set of scalars used for the vector space becomes even more obvious. If we define an endomorphism on a vector space whose scalars are the rational numbers, should we “allow” irrational or complex eigenvalues? You will now recognize our use of the complex numbers in the text for the gross convenience that it is. CBM Change-of-Basis Matrix To create a change-of-basis matrix, it is enough to construct an identity linear transformation relative to a domain and codomain with the specified user bases, which is simply a straight application of Definition CBM. Here we go with two arbitrary bases. b0 = vector(QQ, [-5, 8, 0, 4]) b1 = vector(QQ, [-3, 9, -2, 4]) b2 = vector(QQ, [-1, 4, -1, 2]) b3 = vector(QQ, [-1, 2, 0, 1]) B = [b0, b1, b2, b3] Chapter R 126 U = (QQ^4).subspace_with_basis(B) c0 = vector(QQ, [ 0, 2, -7, 5]) c1 = vector(QQ, [-1, 2, -1, 4]) c2 = vector(QQ, [ 1, -3, 5, -7]) c3 = vector(QQ, [ 1, 1, -8, 3]) C = [c0, c1, c2, c3] V = (QQ^4).subspace_with_basis(C) x1, x2, x3, x4 = var(’x1, x2, x3, x4’) id_symbolic(x1, x2, x3, x4) = [x1, x2, x3, x4] S = linear_transformation(U, V, id_symbolic) CB = S.matrix(side=’right’) CB [ 36 25 8 7] [ 27 34 15 7] [ 35 35 14 8] [-13 -4 0 -2] S.is_invertible() True We can demonstrate that CB is indeed the change-of-basis matrix from B to C, converting vector representations relative to B into vector representations relative to C. We choose an arbitrary vector, x, to experiment with (you could experiment with other possibilities). We use the Sage conveniences to create vector representations relative to the two bases, and then verify Theorem CB. Recognize that x, u and v are all the same vector. x = vector(QQ, [-45, 62, 171, 85]) u = U.coordinate_vector(x) u (-103, -108, 45, 839) v = V.coordinate_vector(x) v (-175, 95, -43, 93) v == CB*u True We can also verify the construction above by building the change-of-basis matrix directly (i.e., without constructing a linear transformation). cols = [V.coordinate_vector(u) for u in U.basis()] M = column_matrix(cols) M Chapter R 127 [ 36 25 8 7] [ 27 34 15 7] [ 35 35 14 8] [-13 -4 0 -2] MRCB Matrix Representation and Change-of-Basis In Sage MR we built two matrix representations of one linear transformation, rel- ative to two different pairs of bases. We now understand how these two matrix representations are related — Theorem MRCB gives the precise relationship with change-of-basis matrices, one converting vector representations on the domain, the other converting vector representations on the codomain. Here is the demonstra- tion. We use MT as the prefix of names for matrix representations, CB as the prefix for change-of-basis matrices, and numerals to distinguish the two domain-codomain pairs. x1, x2, x3, x4 = var(’x1, x2, x3, x4’) outputs = [3*x1 + 7*x2 + x3 - 4*x4, 2*x1 + 5*x2 + x3 - 3*x4, -x1 - 2*x2 + x4] T_symbolic(x1, x2, x3, x4) = outputs U = QQ^4 V = QQ^3 b0 = vector(QQ, [ 1, 1, -1, 0]) b1 = vector(QQ, [-1, 0, -2, 7]) b2 = vector(QQ, [ 0, 1, -2, 4]) b3 = vector(QQ, [-2, 0, -1, 6]) B = [b0, b1, b2, b3] c0 = vector(QQ, [ 1, 6, -6]) c1 = vector(QQ, [ 0, 1, -1]) c2 = vector(QQ, [-2, -3, 4]) C = [c0, c1, c2] d0 = vector(QQ, [ 1, -3, 2, -1]) d1 = vector(QQ, [ 0, 1, 0, 1]) d2 = vector(QQ, [-1, 2, -1, -1]) d3 = vector(QQ, [ 2, -8, 4, -3]) D = [d0, d1, d2, d3] e0 = vector(QQ, [ 0, 1, -3]) e1 = vector(QQ, [-1, 2, -1]) e2 = vector(QQ, [ 2, -4, 3]) E = [e0, e1, e2] U1 = U.subspace_with_basis(B) V1 = V.subspace_with_basis(C) T1 = linear_transformation(U1, V1, T_symbolic) MTBC =T1.matrix(side=’right’) MTBC [ 15 -67 -25 -61] [-75 326 120 298] [ 3 -17 -7 -15] U2 = U.subspace_with_basis(D) V2 = V.subspace_with_basis(E) T2 = linear_transformation(U2, V2, T_symbolic) MTDE = T2.matrix(side=’right’) MTDE Chapter R 128 [ -32 8 38 -91] [-148 37 178 -422] [ -80 20 96 -228] This is as far as we could go back in Section MR. These two matrices represent the same linear transformation (namely T_symbolic), but the question now is “how are these representations related?” We need two change-of-basis matrices. Notice that with different dimensions for the domain and codomain, we get square matrices of different sizes. identity4(x1, x2, x3, x4) = [x1, x2, x3, x4] CU = linear_transformation(U2, U1, identity4) CBDB = CU.matrix(side=’right’) CBDB [ 6 7 -8 1] [ 5 1 -5 9] [-9 -6 10 -9] [ 0 3 -1 -5] identity3(x1, x2, x3) = [x1, x2, x3] CV = linear_transformation(V1, V2, identity3) CBCE = CV.matrix(side=’right’) CBCE [ 8 1 -7] [ 33 4 -28] [ 17 2 -15] Finally, here is Theorem MRCB, relating the the two matrix representations via the change-of-basis matrices. MTDE == CBCE * MTBC * CBDB True We can walk through this theorem just a bit more carefully, step-by-step. We will compute three matrix-vector products, using three vector representations, to demonstrate the equality above. To prepare, we choose the vector x arbitrarily, and we compute its value when evaluted by T_symbolic, and then verify the vector and matrix representations relative to D and E. T_symbolic(34, -61, 55, 18) (-342, -236, 106) x = vector(QQ, [34, -61, 55, 18]) u_D = U2.coordinate_vector(x) u_D (25, 24, -13, -2) Chapter R 129 v_E = V2.coordinate_vector(vector(QQ, [-342, -236, 106])) v_E (-920, -4282, -2312) v_E == MTDE*u_D True So far this is not really new, we have just verified the representation MTDE in the case of one input vector (x), but now we will use the alternate version of this matrix representation, CBCE * MTBC * CBDB, in steps. First, convert the input vector from a representation relative to D to a represen- tation relative to B. u_B = CBDB*u_D u_B (420, 196, -481, 95) Now apply the matrix representation, which expects “input” coordinatized rel- ative to B and produces “output” coordinatized relative to C. v_C = MTBC*u_B v_C (-602, 2986, -130) Now convert the output vector from a representation relative to C to a represen- tation relative to E. v_E = CBCE*v_C v_E (-920, -4282, -2312) It is no surprise that this version of v_E equals the previous one, since we have checked the equality of the matrices earlier. But it may be instructive to see the input converted by change-of-basis matrices before and after being hit by the linear transformation (as a matrix representation). Now we will perform another example, but this time using Sage endomorphisms, linear transformations with equal bases for the domain and codomain. This will allow us to illustrate Theorem SCB. Just for fun, we will do something large. Notice the labor-saving device for manufacturing many symbolic variables at once. [var(’x{0}’.format(i)) for i in range(1, 12)] [x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11] x = vector(SR, [x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11]) A = matrix(QQ, 11, [ [ 146, -225, -10, 212, 419, -123, -73, 3, 219, -100, -57], [ -24, 32, 1, -33, -66, 13, 16, 1, -33, 18, 3], [ 79, -131, -15, 124, 235, -74, -33, -3, 128, -57, -29], [ -1, 13, -16, -1, -27, 3, -5, -4, -9, 6, 2], [-104, 170, 20, -162, -307, 95, 45, 3, -167, 75, 37], [ -16, 59, -19, -34, -103, 27, -10, -1, -51, 27, 8], Chapter R 130 [ 36, -41, -7, 46, 80, -25, -26, 2, 42, -18, -16], [ -5, 0, 1, -4, -3, 2, 6, -1, 0, -2, 3], [ 105, -176, -28, 168, 310, -103, -41, -4, 172, -73, -40], [ 1, 7, 0, -3, -9, 5, -6, -2, -7, 3, 2], [ 74, -141, 4, 118, 255, -72, -23, -1, 133, -63, -26] ]) out = (A*x).list() T_symbolic(x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11) = out V1 = QQ^11 C = V1.basis() T1 = linear_transformation(V1, V1, T_symbolic) MC = T1.matrix(side=’right’) MC [ 146 -225 -10 212 419 -123 -73 3 219 -100 -57] [ -24 32 1 -33 -66 13 16 1 -33 18 3] [ 79 -131 -15 124 235 -74 -33 -3 128 -57 -29] [ -1 13 -16 -1 -27 3 -5 -4 -9 6 2] [-104 170 20 -162 -307 95 45 3 -167 75 37] [ -16 59 -19 -34 -103 27 -10 -1 -51 27 8] [ 36 -41 -7 46 80 -25 -26 2 42 -18 -16] [ -5 0 1 -4 -3 2 6 -1 0 -2 3] [ 105 -176 -28 168 310 -103 -41 -4 172 -73 -40] [ 1 7 0 -3 -9 5 -6 -2 -7 3 2] [ 74 -141 4 118 255 -72 -23 -1 133 -63 -26] Not very interesting, and perhaps even transparent, with a definiton from a matrix and with the standard basis attached to V1 == QQ^11. Let us use a different basis to obtain a more interesting representation. We will input the basis compactly as the columns of a nonsingular matrix. D = matrix(QQ, 11, [[ 1, 2, -1, -2, 4, 2, 2, -2, 4, 4, 8], [ 0, 1, 0, 2, -2, 1, 1, -1, -7, 5, 3], [ 1, 0, 0, -2, 3, 0, -1, -1, 6, -1, -1], [ 0, -1, 1, -1, 3, -2, -3, 0, 5, -8, 2], [-1, 0, 0, 3, -4, 0, 1, 1, -8, 1, 2], [-1, -1, 1, 0, 3, -3, -4, -1, 0, -7, 3], [ 0, 1, 0, 0, 2, 0, 0, -1, 0, -1, 8], [ 0, 0, 0, -1, 0, 0, 0, 1, 5, -4, 1], [ 1, 0, 0, -2, 3, 0, -2, -3, 3, 3, -4], [ 0, -1, 0, 0, 1, -1, -2, -1, 2, -4, 0], [ 1, 0, -1, -2, 0, 2, 2, 0, 5, 3, -1]]) E = D.columns() V2 = (QQ^11).subspace_with_basis(E) T2 = linear_transformation(V2, V2, T_symbolic) MB = T2.matrix(side=’right’) MB [ 2 1 0 0 0 0 0 0 0 0 0] [ 0 2 1 0 0 0 0 0 0 0 0] [ 0 0 2 0 0 0 0 0 0 0 0] [ 0 0 0 2 1 0 0 0 0 0 0] [ 0 0 0 0 2 0 0 0 0 0 0] [ 0 0 0 0 0 -1 1 0 0 0 0] [ 0 0 0 0 0 0 -1 1 0 0 0] [ 0 0 0 0 0 0 0 -1 1 0 0] [ 0 0 0 0 0 0 0 0 -1 1 0] Chapter R 131 [ 0 0 0 0 0 0 0 0 0 -1 1] [ 0 0 0 0 0 0 0 0 0 0 -1] Well, now that is interesting! What a nice representation. Of course, it is all due to the choice of the basis (which we have not explained). To explain the relationship between the two matrix representations, we need a change-of-basis-matrix, and its inverse. Theorem SCB says we need the matrix that converts vector representations relative to B into vector representations relative to C. out = [x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11] id11(x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11) = out L = linear_transformation(V2, V1, id11) CB = L.matrix(side=’right’) CB [ 1 2 -1 -2 4 2 2 -2 4 4 8] [ 0 1 0 2 -2 1 1 -1 -7 5 3] [ 1 0 0 -2 3 0 -1 -1 6 -1 -1] [ 0 -1 1 -1 3 -2 -3 0 5 -8 2] [-1 0 0 3 -4 0 1 1 -8 1 2] [-1 -1 1 0 3 -3 -4 -1 0 -7 3] [ 0 1 0 0 2 0 0 -1 0 -1 8] [ 0 0 0 -1 0 0 0 1 5 -4 1] [ 1 0 0 -2 3 0 -2 -3 3 3 -4] [ 0 -1 0 0 1 -1 -2 -1 2 -4 0] [ 1 0 -1 -2 0 2 2 0 5 3 -1] OK, all set. CB^-1*MC*CB [ 2 1 0 0 0 0 0 0 0 0 0] [ 0 2 1 0 0 0 0 0 0 0 0] [ 0 0 2 0 0 0 0 0 0 0 0] [ 0 0 0 2 1 0 0 0 0 0 0] [ 0 0 0 0 2 0 0 0 0 0 0] [ 0 0 0 0 0 -1 1 0 0 0 0] [ 0 0 0 0 0 0 -1 1 0 0 0] [ 0 0 0 0 0 0 0 -1 1 0 0] [ 0 0 0 0 0 0 0 0 -1 1 0] [ 0 0 0 0 0 0 0 0 0 -1 1] [ 0 0 0 0 0 0 0 0 0 0 -1] Which is MB. So the conversion from a “messy” matrix representation relative to a standard basis to a “clean” representation relative to some other basis is just a similarity transformation by a change-of-basis matrix. Oh, I almost forgot. Where did that basis come from? Hint: find a description of “Jordan Canonical Form”, perhaps in our Second Course in Linear Algebra. CELT Designing Matrix Representations How do we find the eigenvectors of a linear transformation? How do we find pleas- ing (or computationally simple) matrix representations of linear transformations? Theorem EER and Theorem SCB applied in the context of Theorem DC can answer both questions. Here is an example. x1, x2, x3, x4, x5, x6 = var(’x1, x2, x3, x4,x5,x6’) outputs = [ 9*x1 - 15*x2 - 7*x3 + 15*x4 - 36*x5 - 53*x6, 24*x1 - 20*x2 - 9*x3 + 18*x4 - 24*x5 - 78*x6, Chapter R 132 8*x1 - 6*x2 - 3*x3 + 6*x4 - 6*x5 - 26*x6, -12*x1 - 9*x2 - 3*x3 + 13*x4 - 54*x5 - 24*x6, -8*x1 + 6*x2 + 3*x3 - 6*x4 + 6*x5 + 26*x6, -4*x1 - 3*x2 - x3 + 3*x4 - 18*x5 - 4*x6] T_symbolic(x1, x2, x3, x4, x5, x6) = outputs T1 = linear_transformation(QQ^6, QQ^6, T_symbolic) M1 = T1.matrix(side=’right’) M1 [ 9 -15 -7 15 -36 -53] [ 24 -20 -9 18 -24 -78] [ 8 -6 -3 6 -6 -26] [-12 -9 -3 13 -54 -24] [ -8 6 3 -6 6 26] [ -4 -3 -1 3 -18 -4] Now we compute the eigenvalues and eigenvectors of M1. Since M1 is diagonaliz- able, we can find a basis of eigenvectors for use as the basis for a new representation. ev = M1.eigenvectors_right() ev [(4, [ (1, 6/5, 2/5, 4/5, -2/5, 1/5) ], 1), (0, [ (1, 9/7, 4/7, 3/7, -3/7, 1/7) ], 1), (-2, [ (1, 7/5, 2/5, 3/5, -2/5, 1/5) ], 1), (-3, [ (1, 3, 1, -3/2, -1, -1/2) ], 1), (1, [ (1, 0, 0, 3, 0, 1), (0, 1, 1/3, -2, -1/3, -2/3) ], 2)] evalues, evectors = M1.eigenmatrix_right() B = evectors.columns() V = (QQ^6).subspace_with_basis(B) T2 = linear_transformation(V, V, T_symbolic) M2 = T2.matrix(’right’) M2 [ 4 0 0 0 0 0] [ 0 0 0 0 0 0] [ 0 0 -2 0 0 0] [ 0 0 0 -3 0 0] [ 0 0 0 0 1 0] [ 0 0 0 0 0 1] The eigenvectors that are the basis elements in B are the eigenvectors of the lin- ear transformation, relative to the standard basis. For different representations the eigenvectors take different forms, relative to other bases. What are the eigenvectors of the matrix representation M2? Notice that the eigenvalues of the linear transformation are totally independent of the representation. So in a sense, they are an inherent property of the linear transformation. Chapter R 133 You should be able to use these techniques with linear transformations on ab- stract vector spaces — just use a mental linear transformation transforming the abstract vector space back-and-forth between a vector space of column vectors of the right size. SUTH4 Sage Under The Hood, Round 4 We finally have enough theorems to understand how Sage creates and manages linear transformations. With a choice of bases for the domain and codomain, a linear transformation can be represented by a matrix. Every interesting property of the linear transformation can be computed from the matrix representation, and we can convert between representations (of vectors and linear transformations) with change-of-basis matrices, similarity and matrix multiplication. So we can understand the theory of linear algebra better by experimenting with the assistance of Sage, and the theory of linear algebra helps us understand how Sage is designed and functions. A virtuous cycle, if there ever was one. Keep it going.