Linear Algebra

This page supports the main stream of the web site by providing/reminding some information regading Linear Algebra. We demonstrate capabilities of MuPad for this topic.

LU-factorization¶

Recall that in Gaussian Elimination, row operations are used to change the coefficient matrix to an upper triangular matrix. The solution is then found by back substitution, starting from the last equation in the reduced system. In Gauss-Jordan Reduction, row operations are used to diagonalize the coefficient matrix, and the answers are read directly.

The goal of this section is to identify Gaussian elimination with LU factorization. The original m-by-n matrix A becomes the product of two or more special matrices that are actually triangular matrices. Namely, the factorization comes from elimination: \( {\bf A} = {\bf L}\,{\bf U} , \) where L is lower triangular m-by-m matrix and U is upper m-by-n triangular matrix. Such representation is called and LU-decomposition of LU-factorization. Computers usually solve square systems of linear equations using the LU decomposition, and it is also a key step when inverting a matrix, or computing the determinant of a matrix. The LU decomposition was introduced by a Polish astronomer, mathematician, and geodesist Tadeusz Banachiewicz (1882--1954) in 1938.

For a single linear system \( {\bf A}\,{\bf x} = {\bf b} \) of n equations in n unknowns, the methods LU-decomposition and Gauss-Jordan elimination differ in bookkeeping but otherwise involve the same number of flops. However, LU-factorization has the following advantages:

• Gaussian elimination and Gauss--Jordan elimination both use the augmented matrix \( \left[ {\bf A} \, | \, {\bf b} \right] ,\) so b must be known. In contrast, LU-decomposition uses only matrix A, so once that factorization is complete, it can be applied to any vector b.
• For large linear systems in which computer memory is at a premium, one can dispense with the storage of the 1's and zeroes that appear on or below the main diagonal of U, since those entries are known. The space that this opens up can then be used to store the entries of L.

Not every square matrix has an LU-factorization. However, if it is possible to reduce a square matrix A to row echelon form by Gaussian elimination without performing any row interchanges, then A will have an LU-decomposition, though it may not be unique.

LU-decomposition:
Step 1: rewrite the system of algebraic equations \( {\bf A} {\bf x} = {\bf b} \) as

Step 2: Define a new \( n\times 1 \) matrix y (which is actually a vector) by Step 3: Rewrite the given equation as \( {\bf L} {\bf y} = {\bf b} \) and solve this sytem for y.
Step 4: Substitute y into the equation \( {\bf U} {\bf x} = {\bf y} \) and solve for x.

Procedure for constructing LU-decomposition:
Step 1: Reduce \( n \times n \) matrix A to a row echelon form U by Gaussian elimination without row interchanges, keeping track of the multipliers used to introduce the leading coefficients (usually 1's) and multipliers used to introduce the zeroes below the leading coefficients.
Step 2: In each position along the main diagonal L, place the reciprocal of the multiplier that introduced the leading 1 in that position of U.
Step 3: In each position below the main diagonal of L, place the negative of the multiplier used to introduce the zero in that position of U.
Step 4: For the decomposition \( {\bf A} = {\bf L} {\bf U} . \)

Recall the elementary operations on the rows of a matrix that are equivalent to premultiplying by an elementary matrix E:

(1) multiplying row i by a nonzero scalar α , denoted by \( {\bf E}_i (\alpha ) ;\)
(2) adding β times row j to row i, denoted by \( {\bf E}_{ij} (\beta ) \) (here β is any scalar), and
(3) interchanging rows i and j, denoted by \( {\bf E}_{ij} \) (here \( i \ne j \) ), called elementary row operations of types 1,2 and 3, respectively.

Correspondingly, a square matrix E is called an elementary matrix if it can be obtained from an identity matrix by performing a single elementary operation.

Theorem: Every elementary matrix is invertible, and the inverse is also an elementary matrix. ■

Example. Illustrations for m = 4:

Correspondingly, we have their inverses: Multiplying by an arbitrary matrix we obtain

Example. We find LU-factorization of

To achieve this, we reduce A to a row echelon form U using Gaussian elimination and then calculate L by inverting the product of elementary matrices.

Reduction	Row operation	Elementary matrix	Inverse matrix
\( \begin{bmatrix} 2&6&2 \\ -3&-8&0 \\ 4&9&2 \end{bmatrix} \)
Step 1:	\( \frac{1}{2} \times \,\mbox{row 1} \)	\( {\bf E}_1 = \begin{bmatrix} 1/2&0&0 \\ 0&1&0 \\ 0&0&1 \end{bmatrix} \)	\( {\bf E}_1^{-1} = \begin{bmatrix} 2&0&0 \\ 0&1&0 \\ 0&0&1 \end{bmatrix} \)
\( \begin{bmatrix} 1&3&1 \\ -3&-8&0 \\ 4&9&2 \end{bmatrix} \)
Step 2	\( 3 \times \,\mbox{row 1} + \mbox{row 2}\)	\( {\bf E}_2 = \begin{bmatrix} 1&0&0 \\ 3&1&0 \\ 0&0&1 \end{bmatrix} \)	\( {\bf E}_2^{-1} = \begin{bmatrix} 1&0&0 \\ -3&1&0 \\ 0&0&1 \end{bmatrix} \)
\( \begin{bmatrix} 1&3&1 \\ 0&1&3 \\ 4&9&2 \end{bmatrix} \)
Step 3	\( -4 \times \,\mbox{row 1} + \mbox{row 3}\)	\( {\bf E}_3 = \begin{bmatrix} 1&0&0 \\ 0&1&0 \\ -4&0&1 \end{bmatrix} \)	\( {\bf E}_3^{-1} = \begin{bmatrix} 1&0&0 \\ 0&1&0 \\ 4&0&1 \end{bmatrix} \)
\( \begin{bmatrix} 1&3&1 \\ 0&1&3 \\ 0&-3&-2 \end{bmatrix} \)
Step 4	\( 3 \times \,\mbox{row 2} + \mbox{row 3}\)	\( {\bf E}_4 = \begin{bmatrix} 1&0&0 \\ 0&1&0 \\ 0&3&1 \end{bmatrix} \)	\( {\bf E}_4^{-1} = \begin{bmatrix} 1&0&0 \\ 0&1&0 \\ 0&-3&1 \end{bmatrix} \)
\( \begin{bmatrix} 1&3&1 \\ 0&1&3 \\ 0&0&7 \end{bmatrix} \)
Step 5	\( \frac{1}{7} \times \,\mbox{row 3} \)	\( {\bf E}_5 = \begin{bmatrix} 1&0&0 \\ 0&1&0 \\ 0&0&1/7 \end{bmatrix} \)	\( {\bf E}_5^{-1} = \begin{bmatrix} 1&0&0 \\ 0&1&0 \\ 0&0&7 \end{bmatrix} \)
\( \begin{bmatrix} 1&3&1 \\ 0&1&3 \\ 0&0&1 \end{bmatrix} = {\bf U}\)

Now we find the lower triangular matrix: so

Example. Consider the \( 3\times 4 \) matrix

and consider the elementary matrix which results from adding \( -1 \) times the first row of \( {\bf I}_3 \) to the thrid row. The product \( {\bf E}_{31} (-1)\, {\bf A} \) is which is precisely the matrix that results when we add \( -1 \) times the first row of A to the third row. Next, we multiply by matrix \( {\bf E}_{21} (-2) = \begin{bmatrix} 1&0&0 \\ -2&1&0 \\ 0&0&1 \end{bmatrix} . \) This yields Finally, we multiply by \( {\bf E}_{32} (4) = \begin{bmatrix} 1&0&0 \\ 0&1&0 \\ 0&4&1 \end{bmatrix} \) to obtain Hence \( {\bf A} = {\bf E}_{31}^{-1} \,{\bf E}_{21}^{-1} \, {\bf E}_{32}^{-1} \, {\bf U} . \) Evaluation the product of three matrices we obtain the lower triangular matrix L, so \( {\bf A} = {\bf L} \,{\bf U} \) is the required LU-factorization. ■

PLU Factorization

So far, we tried to represent a square nonsingular matrix A as a product of a lower-triangular matrix L and an upper triangular matrix U: \( {\bf A} = {\bf L}\,{\bf U} .\) When this is possible we say that A has an LU-decomposition (or factorization). It turns out that this factorization (when it exists) is not unique. If L has 1's on it's diagonal, then it is called a Doolittle factorization. If U has 1's on its diagonal, then it is called a Crout factorization. When L is an unitary matrix, it is called a Cholesky decomposition. Moreover, Gaussian elimination in its pure form may be unstable.

While LU-decomposition is a useful computational tool, but this does not work for every matrix. Consider even the simple example with matrix

Another cause of instability in Gaussian elimination is when encountering relatively small pivots. Gaussian elimination is guaranteed to succeed if row or column interchanges are used in order to avoid zero pivots when using exact calculations. But, when computing the LU factorization numerically, this is not necessarily true. When a calculation is undefined because of division by zero, the same calculation will suffer numerical difficulties when there is division by a nonzero number that is relatively small. When computing the factors L and U, the process does not fail in this case because there is no division by zero. When these computations are performed in floating point arithmetic, the number \( 3 - 2\cdot 10^{22} \) is not represented exactly but will be rounded to the nearest floating point number which we will say is \( 2\cdot 10^{22} . \) This means that our matrices are now floating point matrices L' and U': When we compute L'U', we get The answer should be equal to A, but obviously that is not the case. The 3 in position (2,2) of matrix A is now 0. Also, when trying to solve a system such as \( {\bf A} \, {\bf x} = {\bf b} \) using the LU factorization, the factors L'U' would not give you a correct answer. The LU factorization was a stable computation but not backward stable.

Recall that an \( n \times n \) permutation matrix P is a matrix with precisely one entry whose value is "1" in each column and row, and all of whose other entries are "0." The rows of P are a permutation of the rows of the identity matrix. Since not every matrix has LU decomposition, we try to find a permulation matrix P such that PA has LU factorization: \( {\bf P}\,{\bf A} = {\bf L}\,{\bf U} , \) where L and U are again lower and upper triangular matrices, and P is a permutation matrix. In this case, we say that A has a PLU factorization. It is also referred to as the LU factorization with Partial Pivoting (LUP) with row permutations only. An LU factorization with full pivoting involves both row and column permutations, \( {\bf P}\,{\bf A}\, {\bf Q} = {\bf L}\,{\bf U} , \) where L and U, and P are defined as before, and Q is a permutation matrix that reorders the columns of A.

For a permutation matrix P, the product PA is a new matrix whose rows consists of the rows of A rearranged in the new order. Note that a product of permutation matrices is a permutation matrix. Every permutation matrix is an orthogonal matrix: \( {\bf P}^{-1} = {\bf P}^{\mathrm T} . \)

Example. Let P be a permutation matrix that interchange rows 1 and 2 and also interchange rows 3 and 4:

However,

Example. Consider a 3-by-4 matrix \( {\bf A} = \begin{bmatrix} 1&2&-3&1 \\ 2&4&0&7 \\ -1&3&2&0 \end{bmatrix} \) for which pure LU-factorization does not exist. Multiplying A by the product of two matrices

we obtain To interchange last two rows, we multiply by the matrix \( {\bf P} = \begin{bmatrix} 1&0&0 \\ 0&0&1 \\ 0&1&0 \end{bmatrix} . \) This yields From this equation, we find the inverse matrix which is not a triangular matrix. Therefore, we get PLU decomposition: Note that \( {\bf P}^2 = {\bf I} , \) the identity matrix, and \( {\bf P}\, {\bf A} = {\bf L}\,{\bf U} . \)

Theorem: If A is a nonsingular matrix, then there exists a permutation matrix P so that PA has an LU-factorization. ■

Pivoting for LU factorization is the process of systematically selecting pivots for Gaussian elimination during the LU-decomposition of a matrix. The LU factorization is closely related to Gaussian elimination, which is unstable in its pure form. To guarantee the elimination process goes to completion, we must ensure that there is a nonzero pivot at every step of the elimination process. This is the reason we need pivoting when computing LU decompositions. But we can do more with pivoting than just making sure Gaussian elimination completes. We can reduce roundoff errors during computation and make our algorithm backward stable by implementing the right pivoting strategy. Depending on the matrix A, some LU decompositions can become numerically unstable if relatively small pivots are used. Relatively small pivots cause instability because they operate very similar to zeroes during Gaussian elimination. Through the process of pivoting, we can greatly reduce this instability by ensuring that we use relatively large entries as our pivot elements. This prevents large factors from appearing in the computed L and U, which reduces roundoff errors during computation.

The goal of partial pivoting is to use a permutation matrix to place the largest entry of the first column of the matrix at the top of that first column. For an \( n \times n \) matrix A, we scan n rows of the first column for the largest value. At step k of the elimination, the pivot we choose is the largest of the n- (k+1) subdiagonal entries of column k, which costs O(nk) operations for each step of the elimination. So for an n-by-n matrix, there is a total of \( O(n^2 ) \) comparisons. Once located, this entry is then moved into the pivot position \( A_{k,k} \) on the diagonal of the matrix.

Example. Use PLU factorization to solve the linear system

We want to use a permutation matrix to place the largest entry of the first column into the position (1,1). For this matrix, this means we want 4 to be the first entry of the first column. So we use the permutation matrix Then we use the matrix M to eliminate the bottom two entries of the first column: For latter matrix we don't need pivoting because the middle of the second row is the largest. So we just multiply by Since the inverse of \( {\bf M}_2 {\bf M}_1 \) is the required lower triangular matrix we have then Thus, the solution becomes

Algorithm for PLU-factorization: is based on the following formula for upper triangular matrix \( {\bf U} = \left[ u_{i,j} \right] : \)

with similar formula for \( {\bf L} = \left[ l_{i,j} \right] . \) To ensure that the algorithm is numerically stable when \( u_{j,j} \ll 0, \) a pivoting matrix P is used to re-order A so that the largest element of each column of A gets shifted to the diagonal of A. The formula for elements of L follows:

Step 1: If \( n=1 , \) return \( {\bf P} =1, \quad {\bf L} = 1, \quad {\bf U} = a_{1,2} . \)

Step 2: Find a non-zero element \( a_{i,1} \) in the first row of A. (If not, then the system is not solvable.)

Step 3: Define a permutation matrix \( {\bf Q}_{n \times n} \) such that

• \( q_{i,1} = q_{1,i} =1 . \)
• \( q_{j,j} =1 \) for \( j \ne 1, \ j\ne i . \)
• All the other elements of Q are 0. ■

Step 4: Let \( {\bf B} = {\bf Q} {\bf A} . \) Note that \( b_{1,1} = a_{i,1} \ne 0. \)

Step 5: Write B as \( {\bf B} = \begin{bmatrix} b_{1,1} & {\bf r}^{\mathrm T} \\ {\bf c} & {\bf C}_1 \end{bmatrix} . \)

Step 6: Define \( {\bf B}_1 = {\bf C}_1 - {\bf c} \cdot {\bf r}^{\mathrm T} / b_{1,1} .\) Recursively compute \( {\bf B}_1 = {\bf P}_1 {\bf L}_1 {\bf U}_1 . \)

Step 7:

When employing the complete pivoting strategy we scan for the largest value in the entire submatrix \( {\bf A}_{k:n,k:n} , \) where k is the number of the elimination, instead of just the next subcolumn. This requires \( O \left( (n-k)^2 \right) \) comparisons for every step in the elimination to find the pivot. The total cost becomes \( O \left( n^3 \right) \) comparisons for an n-by-n matrix. For an \( n \times n \) matrix A, the first step is to scan n rows and n columns for the largest value. Once located, this entry is then permuted into the next diagonal pivot position of the matrix. So in the first step the entry is permuted into the (1,1) position of matrix A. We interchange rows exactly as we did in partial pivoting, by multiplying A on the left with a permutation matrix P. Then to interchange columns, we multiply A on the right by another permutation matrix Q. The matrix product PAQ interchanges rows and columns accordingly so that the largest entry in the matrix is in the (1,1) position of A. With complete pivoting, the general equation for L is the same as for partial pivoting, but the equation for U is slightly different. Complete pivoting is theoretically the most stable strategy as it can be used even when a matrix is singular, but it is much slower than partial pivoting.

Example. Consider the matrix A:

In this case, we want to use two different permutation matrices to place the largest entry of the entire matrix into the position \( A_{1,1} .\) For the given matrix, this means we want 9 to be the first entry of the first column. To achieve this we multiply A on the left by the permutation matrix \( {\bf P}_1 \) and multiply on the right by the permutation matrix \( {\bf Q}_1 .\) Then in computing the LU factorization, the matrix M is used to eliminate the bottom two entries of the first column: Then, the permutation matrices P and Q are used to move the largest entry of the submatrix \( {\bf A}_{2:3;2:3} \) into the position \( A_{2,2} : \) As we did before, we are going to use matrix multiplication to eliminate entries at the bottom of the column 2: This upper triangular matrix product is the U matrix. The L matrix comes from the inverse of the product of M and P matrices: To check, we can show that \( {\bf P} {\bf A} {\bf Q} = {\bf L} {\bf U} , \) where \( {\bf P} = {\bf P}_2 {\bf P}_1 \) and \( {\bf Q} = {\bf Q}_1 {\bf Q}_2 . \) ■

Rook pivoting is a third alternative pivoting strategy. For this strategy, in the kth step of the elimination, where \( 1 \le k \le n-1 \) for an n-by-n matrix A, we scan the submarix \( {\bf A}_{k:n;k:n} \) for values that are the largest in their respective row and column. Once we have found our potential pivots, we are free to choose any that we please. This is because the identified potential pivots are large enough so that there will be a very minimal difference in the end factors, no matter which pivot is chosen. When a computer uses rook pivoting, it searches only for the first entry that is the largest in its respective row and column. The computer chooses the first one it finds because even if there are more, they all will be equally effective as pivots. This allows for some options in our decision making with regards to the pivot we want to use.

Example. Consider the matrix A:

When using the partial pivoting technique, we would identify 8 as our first pivot element. When using complete pivoting, we would use 9 as our first pivot element. When rook pivoting we identify entries that are the largest in their respective rows and columns. So following this logic, we identify 8, 9, and 7 all as potential pivots. In this example we will choose 7 as the first pivot. Since we are choosing 7 for our first pivot element, we multiply A on the right by the permutation matrix Q to move 7 into the position (1,1). Note that we will also multiply on the left of A by P, but P is simply the identity matrix so there will be no effect on A. Then we eliminate the entries in row 2 and 3 of column 1 via the matrix If we were employing partial pivoting, our pivot would be 48/7 and if we were using complete pivoting, our pivot would be 62/7. But using the rook pivoting strategy allows us to choose either 48/7 or 62/7 as our next pivot. In this situation, we do not need to permute the matrix \( {\bf M}_1 {\bf P}_1 {\bf A}\,{\bf Q}_1 \) to continue with our factorization because one of our potential pivots is already in pivot postion. To keep things interesting, we will permute the matrix anyway and choose 62/7 to be our next pivot. To achieve this we wil multiply by the matrix \( {\bf P}_2 \) on the left and by the matrix \( {\bf Q}_2 \) on the right. At this point, we want to eliminate the last entry of the second column to complete the factorization. We achieve this by matrix \( {\bf M}_2 :\) After the elimination the matrix is now upper triangular which we recognize as the factor U. The factor L is equal to the matrix product \( {\bf L} = \left( {\bf M}_2 {\bf P}_2 {\bf M}_1 {\bf P}_1 \right)^{-1} . \) You can use Sage to check that PAQ = LU, where \( {\bf P} = {\bf P}_2 {\bf P}_1 \) and \( {\bf Q} = {\bf Q}_1 {\bf Q}_2 \)