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. For instance, let
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':
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 1:
Let P be a permutation matrix that interchange rows 1 and 2 and
also interchange rows 3 and 4:
Example 2:
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
Note that \( {\bf P}^2 = {\bf I} , \) the identity matrix, and \( {\bf P}\, {\bf A} =
{\bf L}\,{\bf U} . \)
End of Example 2
■
Theorem 1:
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 3:
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
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:
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.
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 .\)
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 . \)
End of Example 4
■
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.
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.
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.
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 \)
A. Edelman. The complete pivoting conjecture for Gaussian elimination is false.
The Mathematica Journal, 2:58–61, 1992.
L. V. Foster. The growth factor and efficiency of Gaussian elimination with rook pivoting.
J. Comput.
Appl. Math.
, 86:177–194, 1997. Corrigendum in J. Comput. Appl. Math., 98:177, 1998.
N. I. M. Gould. On growth in Gaussian elimination with complete pivoting.
SIAM J. Matrix Anal. Appl.
,
12(2):354–361, 1991.
N. J. Higham and D. J. Higham. Large growth factors in Gaussian elimination with pivoting.
SIAM J.
Matrix Anal. Appl.
, 10(2):155–164, Apr. 1989.
S. Toledo. Locality of reference in LU decomposition with partial pivoting.
SIAM J. Matrix Anal. Appl.
,
18(4):1065–1081, 1997.