Lecture #07

Scientific Computing Using Python - PHYS:4905 - Fall 2018

Lecture #07 - 9/11/2018 - Prof. Kaaret

These notes borrow from Linear Algebra by Cherney, Denton, Thomas, and Waldron.

Elementary row operations as matrices

Our three elementary row operations can be translated matrices. Let's start with 'Row swap' or exchanging two rows. The following matrix will swap the first and second rows of an augmented matrix.

$\left(\begin{array}{ccc} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \\ \end{array}\right)$

What matrix would swap the first and third rows?

The matrix for scalar multiplication, multiplication of a row by a non-zero constant, is particularly simple. This matrix will multiply the third row by the scalar c and leave the other rows unchanged. In general, one starts with the identity matrix and and changes one element on the diagonal to be not equal to one.

$\left(\begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & c \\ \end{array}\right)$

For row addition, addition of a multiple of one row to another row, one again starts with the identity matrix, but puts in an off diagonal element. What will this matrix do?

$\left(\begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & c \\ 0 & 0 & 1 \\ \end{array}\right)$

The inverse of a matrix

If we want to solve an equation like 7x = 3, we multiply both sides of the equation by 1/7, so

7x = 3
(1/7) 7x = (1/7) 3
x = 3/7
The value 1/7 is the inverse of 7.

Can we do this with matrices?

If M and N are matrices and I is the identity matrix, then we say that N is the inverse of M if NM = I.

Sometimes we write N = M^-1, in analogy with 1/7 = 7^-1.

So, how do you find the inverse of a matrix? Can you divide the identity matrix by M?

Let's go back to our system of equations

$\begin{matrix} x + y = 27 \\ 2x - y = 0\\ \end{matrix} \;\;\; ⟺ \;\;\; \begin{pmatrix} 1 & 1 \\ 2 & -1 \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} = \begin{pmatrix} 27 \\ 0 \end{pmatrix}$

We learned Gaussian elimination so that we could solve this set of equations by transforming the augmented matrix

\left(\begin{array}{cc|l} 1 & 1 & | \, 27 \\ 2 & -1 & | \, 0 \\ \end{array}\right) \;\;\; ⟺ \;\;\; \begin{pmatrix} 1 & 1 \\ 2 & -1 \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} = \begin{pmatrix} 27 \\ 0 \end{pmatrix}

into

\left( \begin{array}{cc|l} 1 & 0 & | \, 9 \\ 0 & 1 & | \, 18 \\ \end{array}\right) \;\;\; ⟺ \;\;\; \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} = \begin{pmatrix} 9 \\ 18 \end{pmatrix}

We did the Gaussian elimination by applying a set of elementary row operations. From the first part of this lecture, we can now write each of those EROs as a matrix. Let's call them E₁, E₂, .E₃, ...

If we call our original matrix M, so

N = \begin{pmatrix} 1 & 1 \\ 2 & -1 \end{pmatrix}

Then by applying our EROs, we get the identity matrix, so

E_n ... E₂ E₁ M = I

Note the order. You need to multiply M by E₁ first, since it is the first operation applied.

Matrix multiplication is associative, you can change the grouping of the multiplications, but it is not commutative, you cannot, in general, swap the order of the multiplications.

Since matrix multiplication is associative, we can re-write the equation above as

(E_n ... E₂ E₁ ) M = I

Comparing this to equation that defines the inverse of M, NM = I, we see that the product of all of the EROs gives the inverse.

N = M^-1 = E_n ... E₂ E₁

Isn't that nice?

Note that we were able to transform the left part of the augmented matrix to the identity matrix only when the corresponding set of linear equations has a unique set of solutions. This wasn't possible for all matrices. In the same way, not all matrices have an inverse. Indeed, a matrix M is invertible only if its reduced row echelon form is an identity matrix.

Interestingly, you can use the EROs to build up M if you like. All of the EROs matrices that we defined above have inverses. You can define a new ERO to undo the effect of any given ERO, so this must be true. Now let's try

M = E₁^-1 E₂^-1 ... E_n^-1

By the definition of the inverse,

M^-1 M = I

If we plug in for M and M^-1, we find

E_n ... E₂ E₁E₁^-1 E₂^-1 ... E_n^-1 = I

Remember that we can associate the matrix multiplications any way we like. A nice way to associate them is to start at the middle where we have E₁E₁^-1 . This, of course, equals the identity matrix. Multiplication by the identity matrix gives back the original matrix, or we can just cancel out the identity matrix from the equation, just like you remove factors that multiply out to 1 in algebra. The middle of the equation will then be E₂ E₂^-1 . We can repeat this procedure until we have only the identity matrix on the left hand side of the equation, leaving I = I, which is always true, so our original guess about writing M in terms of the inverses of the EROs must also be true.

LU Factorization

In addition to finding inverses, we can also factor a matrix. One form of a matrix that is useful for numerical computation is the LU factorization. The matrix M is written as the produce M = LU where L is a lower triangular matrix, meaning that only the lower elements along or below the diagonal are non-zero, and U is an upper triangular matrix meaning only upper elements along or above the diagonal are non-zero. For example,

$M = LU = \begin{pmatrix} 3 & 0 & 0 \\ 2 & 5 & 0 \\ 1 & 7 & 4 \end{pmatrix} \begin{pmatrix} 2 & 0 & 9 \\ 0 & 1 & 5 \\ 0 & 0 & 6 \end{pmatrix}$

One can also do LDU factorization where the diagonal elements of the L and U matrices are 1 and there is a diagonal matrix in the middle.

To do a factorization, one does Gaussian elimination and then groups the EROs into blocks. The L matrix is the product of the inverses of EROs which zero out below the diagonal by row addition, D is the product of inverses of EROs which set the diagonal elements to 1 by row multiplication, and U is the product of inverses of EROs which zero out above the diagonal by row addition. If row swaps are needed, then another matrix is needed in the factorization. Note that the EROs need to be done in the correct order, since matrix multiplication is not communitative. There are algorithms, which are basically variations of Gaussian elimination, to figure out the EROs.

Sets of solutions to sets of linear equations

A single linear equation can have:
- no solutions, e.g. 0x = 9 has no solutions,
- one solution, e.g. 3x = 9 has one solution x = 3,.
- or infinitely many solutions., i.e. 0x = 0 is solved by any real number.

A pair of linear equations can have also have no solutions, one solution, or infinitely many solutions, but there are two types of infinitely many solutions. One could have the equations x + 2y = 3 and 2x + 4y = 6 for which any real number can be a solution for x as long as the cooresponding y = (3 - x)/2. Or one could have 0x + 0y = 0 and 0x + 0y = 0 in which any pair of real numbers is a solution.

If you think of this geometrically, then the set of possible values of (x, y) are a two dimensional space. The case of one solution is a point in the space, the pair of equations x + 2y = 3 and 2x + 4y = 6 have a set of solutions that form a line in the space, and the zero equations fill the whole space. If we have n variables, then they form an n-dimensional space and the solutions can range from a point in the space, to a line, to a plane, ..., to filling up the whole space.

We can write our system of linear equations in matrix form as Mx = v. If the system of linear equations has multiple solutions, meaning at least one line of solutions, then it must be true that the equation Mx = 0 has some solution x^H such that Mx^H = 0. The line of solutions can then be written in the form x = x^P + µx^H where is a real number. The equation Mx = 0 his called a homogeneous equation and its solutions are called the homogeneous solutions. If the set of a solutions is a plane, then must be two homogenous solutions that one can use to define the plane.

Assignment

The assignment for today's class involves converting EROs to matrix form and using them to find an inverse matrix and also writing a general program to do Gaussian elimination.

HW #7 is due at the beginning of class on Tuesday, September 18.
https://homepage.physics.uiowa.edu/~pkaaret/2018f_p4905/hw07.html