AMTH247 Lecture 9

Linear Equations III

Reading: Heath §2.4.5

Contents

• Pivoting
  • The Need for Pivoting
  • Partial Pivoting
  • Permutation Matrices
  • LU Factorization with Partial Pivoting
  • The LU Factorization Algorithm
  • Scilab
  
  
• About this document ...

Pivoting

The Need for Pivoting

At each stage of the process of LU factorization we have a partially reduced matrix of the form

$$\begin{bmatrix} a_{11} & a_{12} & \ldots & a_{1k} & \ldots & a_{... ...ts & \vdots \\ 0 & 0 & \ldots & a_{nk} & \ldots & a_{nn} \\ \end{bmatrix} ,$$

and we want to set to zero the elements in column $k$ below $a_{kk}$. The method for doing this, multiplying by the elimination matrix

$$\mathbf{M}_k = \begin{bmatrix} 1 & \ldots & 0 & 0 & \ldots & 0... ...dots & \ddots & \vdots \\ 0 & \ldots & -m_{n} & 0 & \ldots & 1 \end{bmatrix}$$

where

$$m_i = -\frac{a_{ik}}{a_{kk}}, \quad i = k+1, \ldots, n$$

fails when the divisor $a_{kk}$ is zero. The solution is obvious, if $a_{kk} = 0$ swap the $k$-th with some other row $i$ with $i > k$, for which $a_{ik} \neq 0$. This is always possible unless $a_{ik} = 0$ for all $i = k, \ldots, n$, in which case the matrix must be singular. Even when $a_{kk}$ is not zero, there may be advantages in swapping rows to reduce rounding error. We saw in the previous lecture that this is what Scilab's LU factorization does.

Partial Pivoting

As has been mentioned, the point of pivoting, i.e. the interchange of rows, is to to minimize the effect of rounding errors¹. There are a number of approaches to pivoting, the most important being partial pivoting. In partial pivoting, when we reach the stage

$$\begin{bmatrix} a_{11} & a_{12} & \ldots & a_{1k} & \ldots & a_{... ...ts & \vdots \\ 0 & 0 & \ldots & a_{nk} & \ldots & a_{nn} \\ \end{bmatrix} ,$$

we choose as the pivot element the element $a_{ik}$, $i \geq k$, of largest magnitude, that is the largest element on or below the diagonal in the $k$th column, and swap row $i$ with row $k$. Intuitively the justification for this is that by choosing the pivot element as large as possible, the multipliers are as small as possible and previous rounding errors will be not be amplified in the multiplication by the elimination matrix.

Example:

We will use the same example as in Lecture 8:

$$\textbf{A} = \begin{bmatrix}1 & 2 & 3 \\ 2 & 3 & 1 \\ 3 & 1 & 2 \end{bmatrix} .$$

The element of largest magnitude in the first column is 3, so we swap the first and third rows

$$\textbf{A}' = \begin{bmatrix}3 & 1 & 2 \\ 2 & 3 & 1 \\ 1 & 2 & 3 \end{bmatrix} .$$

The elimination matrix is

$$\mathbf{M}_1 = \begin{bmatrix}1 & 0 & 0 \\ -\frac{2}{3} & 1 & 0 \\ -\frac{1}{3} & 0 & 1 \end{bmatrix}$$

and

$$\mathbf{A}_1 = \mathbf{M}_1 \mathbf{A}' = \begin{bmatrix}3 & 1 ... ...rac{1}{3} & -\frac{1}{3} \\ 0 & 1 \frac{2}{3} & 2 \frac{1}{3} \end{bmatrix} .$$

The largest element on or below the diagonal in the second column is $2 \frac{1}{3}$ so no row interchange is needed at this step. The elimination matrix is

$$\mathbf{M}_2 = \begin{bmatrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -\frac{5}{7} & 1 \end{bmatrix}$$

and

$$\mathbf{A}_2 = \mathbf{M}_2 \mathbf{A} _1= \begin{bmatrix}3 & 1... ... \\ 0 & 2 \frac{1}{3} & -\frac{1}{3} \\ 0 & 0 & 2 \frac{4}{7} \end{bmatrix} .$$

This is the same upper triangular matrix which was computed by Scilab at the end of Lecture 8.

Permutation Matrices

A permutation matrix is a $n \times n$ matrix $\mathbf{P}$ with exactly one 1 in each row and each column and zeros elsewhere, or, equivalently, an identity matrix with its rows and columns permuted. Some properties of permutation matrices:

1. The inverse of a permutation matrix $\mathbf{P}$ is its transpose $\mathbf{P}^T$.
2. The product of two permutation matrices is again a permutation matrix.
3. If $\mathbf{P}$ is a permutation matrix and $\mathbf{A}$ is any $n \times n$ matrix then
   1. The matrix product $\mathbf{P} \mathbf{A}$ permutes the rows of $\mathbf{A}$.
   2. The matrix product $\mathbf{A} \mathbf{P}$ permutes the columns of $\mathbf{A}$.

Example:

Let $\mathbf{P}$ be the permutation matrix

$$\mathbf{P} = \begin{bmatrix}0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix}$$

and let

$$\mathbf{A} = \begin{bmatrix}1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{bmatrix} .$$

Then

$$\mathbf{P} \mathbf{A} = \begin{bmatrix}7 & 8 & 9 \\ 1 & 2 & 3 \\ 4 & 5 & 6 \end{bmatrix}$$

permutes the rows of $\mathbf{A}$ in the order $(3,1,2)$ and

$$\mathbf{A} \mathbf{P} = \begin{bmatrix}2 & 3 & 1 \\ 5 & 6 & 4 \\ 8 & 9 & 7 \end{bmatrix}$$

permutes the columns of $\mathbf{A}$ in the order $(3,1,2)$.

LU Factorization with Partial Pivoting

When performing LU factorization with partial pivoting we must keep track of the row interchanges we perform. Row interchanges can be performed by multiplying by a permutation matrix. To interchange rows $i$ and $j$ of a matrix $\mathbf{A}$ we multiply by the permutation matrix $\mathbf{P}_{ij}$ obtained by swapping rows $i$ and $j$ of the identity matrix. For example, the following $6 \times 6$ permutation matrix swaps the second and fourth rows:

$$\mathbf{P}_{24} = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 \\ 0 ... ... & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix}$$

These particular permutation matrices are symmetric, i.e. satisfy $$\mathbf{P}^T = \mathbf{P}$$, and so are their own inverse, $$\mathbf{P}^{-1} = \mathbf{P}$$. The effect of partial pivoting is to multiply by a permutation matrix before multiplying by the elimination matrix at each step. So instead of performing the elimination by

$$\mathbf{U} = \mathbf{M}_{n-1} \ldots \mathbf{M}_2 \mathbf{M}_1 \mathbf{A} ,$$

where the $\mathbf{M}_j$ are elimination matrices, we compute instead

$$\mathbf{U} = \mathbf{M}_{n-1} \mathbf{P}_{n-1} \ldots \mathbf{M}_2 \mathbf{P}_2 \mathbf{M}_1 \mathbf{P}_1 \mathbf{A} ,$$

where the $\mathbf{P}_j$ are permutation matrices. Now the inverse of the matrix

$$\mathbf{M} = \mathbf{M}_{n-1} \mathbf{P}_{n-1} \ldots \mathbf{M}_2 \mathbf{P}_2 \mathbf{M}_1 \mathbf{P}_1$$

is

$$\mathbf{L} = \mathbf{P}_1 \mathbf{L}_1 \mathbf{P}_2 \mathbf{L}_2 \ldots \mathbf{P}_{n-1} \mathbf{L}_{n-1}$$

where $\mathbf{L}_j$ is the inverse of $\mathbf{M}_j$ and, as mentioned above, the $\mathbf{P}_j$ are their own inverses. Once more we have the factorization

$$\mathbf{A} = \mathbf{L} \mathbf{U}$$

with $\mathbf{L}$ and $\mathbf{U}$ as above, but while $\mathbf{U}$ is upper triangular, the matrix $\mathbf{L}$ is, in general, no longer lower triangular because of the permutation matrices in its construction.

Example:

Starting with the matrix

$$\textbf{A} = \begin{bmatrix}1 & 2 & 3 \\ 2 & 3 & 1 \\ 3 & 1 & 2 \end{bmatrix} .$$

we computed

$$\mathbf{U} = \begin{bmatrix}3 & 1 & 2 \\ 0 & 2 \frac{1}{3} & -\frac{1}{3} \\ 0 & 0 & 2 \frac{4}{7} \end{bmatrix} .$$

by first swapping the first and third rows, corresponding to the permutation matrix

$$\mathbf{P}_1 = \begin{bmatrix}0 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 & 0 \end{bmatrix} ,$$

then applying the elimination matrix

$$\mathbf{M}_1 = \begin{bmatrix}1 & 0 & 0 \\ -\frac{2}{3} & 1 & 0 \\ -\frac{1}{3} & 0 & 1 \end{bmatrix} .$$

There is now swapping of rows at the second step, so

$$\mathbf{P}_2 = \begin{bmatrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} .$$

The second elimination matrix was

$$\mathbf{M}_2 = \begin{bmatrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -\frac{5}{7} & 1 \end{bmatrix} ,$$

so the matrix $\mathbf{U}$ is

$$\mathbf{U} = \mathbf{M}_2 \mathbf{P}_2 \mathbf{M}_1 \mathbf{P}_1 \mathbf{A} ,$$

The matrix $\mathbf{L}$ is then given by

$$\mathbf{L} = \mathbf{P}_1 \mathbf{L}_1 \mathbf{P}_2 \mathbf{L}_2$$

$$= \begin{bmatrix}0 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 & 0 \end{bmatr... ... \begin{bmatrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & \frac{5}{7} & 1 \end{bmatrix}$$

$$= \begin{bmatrix}\frac{1}{3} & \frac{5}{7} & 1 \\ \frac{2}{3} & 1 & 0 \\ 1 & 0 & 0 \end{bmatrix}$$

This is exactly the result given by Scilab.

The LU Factorization Algorithm

With partial pivoting the matrix $\mathbf{L}$ is not lower triangular, however it is a row permutation of a lower triangular matrix. More precisely we compute a factorization

$$\mathbf{A} = \mathbf{L} \mathbf{U}$$

such that $\mathbf{U}$ is upper triangular and there is a permutation matrix $\mathbf{P}$ such that $$\mathbf{L}' = \mathbf{P} \mathbf{L}$$ is lower triangular. Multiplying the equation above by the permutation matrix $\mathbf{P}$ we have, equivalently,

$$\mathbf{P} \mathbf{A} = \mathbf{P} \mathbf{L} \mathbf{U}$$

$$= \mathbf{L}' \mathbf{U}$$

where $$\mathbf{L}'$$ is lower triangular. The permutation matrix $\mathbf{P}$ is given by

$$\mathbf{P} = \mathbf{P}_{n-1} \ldots \mathbf{P}_2 \mathbf{P}_1$$

where the matrices $\mathbf{P}_j$ are the permutation matrices associated with pivoting introduced in previous section. Suppose we now have a LU factorization

$$\mathbf{L} \mathbf{U} = \mathbf{P} \mathbf{A}$$

and we wish to solve a linear system

$$\mathbf{A} \mathbf{x} = \mathbf{b} .$$

Multiply both sides by $\mathbf{P}$ to give

$$\mathbf{P} \mathbf{A} \mathbf{x} = \mathbf{P} \mathbf{b}$$
or

$$\mathbf{L} \mathbf{U} \mathbf{x} = \mathbf{P} \mathbf{b} .$$

This system can be solved by successive substitution:

$$\mathbf{L} \mathbf{z} = \mathbf{P} \mathbf{b}$$

$$\mathbf{U} \mathbf{x} = \mathbf{z} .$$

Example:

Previously we computed

$$\mathbf{L} = \begin{bmatrix}\frac{1}{3} & \frac{5}{7} & 1 \\ \f... ... \\ 0 & 2 \frac{1}{3} & -\frac{1}{3} \\ 0 & 0 & 2 \frac{4}{7} \end{bmatrix} .$$

and the permutation matrices used were

$$\mathbf{P}_1 = \begin{bmatrix}0 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 &... ...hbf{P}_2 = \begin{bmatrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} ,$$

Now

$$\mathbf{P} = \mathbf{P}_2 \mathbf{P}_1$$

$$= \begin{bmatrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix}0 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 & 0 \end{bmatrix}$$

$$= \begin{bmatrix}0 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 & 0 \end{bmatrix}$$

and the required lower triangular form of $\mathbf{L}$ is

$$\mathbf{L} = \mathbf{P} \mathbf{L}$$

$$= \begin{bmatrix}0 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 & 0 \end{bmatrix... ...rac{1}{3} & \frac{5}{7} & 1 \\ \frac{2}{3} & 1 & 0 \\ 1 & 0 & 0 \end{bmatrix}$$

$$= \begin{bmatrix}1 & 0 & 0 \\ \frac{2}{3} & 1 & 0 \\ \frac{1}{3} & \frac{5}{7} & 1 \end{bmatrix}$$

Thus we have the factorization

$$\mathbf{L} \mathbf{U} = \mathbf{P} \mathbf{A}$$

with

$$\mathbf{L} =\begin{bmatrix}1 & 0 & 0 \\ \frac{2}{3} & 1 & 0 \\ ... ...mathbf{P} = \begin{bmatrix}0 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 & 0 \end{bmatrix}$$

To solve the linear system

$$\mathbf{A} \mathbf{x} = \mathbf{b}$$

we compute

$$\mathbf{b}' = \mathbf{P} \mathbf{b} = \begin{bmatrix}0 & 0 & 1 ... ...{bmatrix}1 \\ 1 \\ 1 \end{bmatrix} = \begin{bmatrix}1 \\ 1 \\ 1 \end{bmatrix}$$

and then

1. Solve $$\mathbf{L} \mathbf{z} = \mathbf{b}'$$ for $\mathbf{z}$. The solution is
   $$\mathbf{z} = \begin{bmatrix}1 \\ 1/3 \\ 3/7 \end{bmatrix}$$

2. Solve $\mathbf{U} \mathbf{x} = \mathbf{z}$ for $\mathbf{x}$. The solution is, as before,
   $$\mathbf{x} = \begin{bmatrix}1/6 \\ 1/6 \\ 1/6 \end{bmatrix}$$

Scilab

There is a second form of the lu command in Scilab:

    [l,u,p] = lu(a)

which computes the LU factorization we have described above. For our example

-->a = [1 2 3
-->2 3 1
-->3 1 2]
 a  =
 
!   1.    2.    3. !
!   2.    3.    1. !
!   3.    1.    2. !
 
-->[l,u,p] = lu(a)
 p  =
 
!   0.    0.    1. !
!   0.    1.    0. !
!   1.    0.    0. !
 u  =
 
!   3.    1.           2.        !
!   0.    2.3333333  - 0.3333333 !
!   0.    0.           2.5714286 !
 l  =
 
!   1.           0.           0. !
!   0.6666667    1.           0. !
!   0.3333333    0.7142857    1. !

which agrees with our calculations.

About this document ...

AMTH247 Lecture 9

Linear Equations III

This document was generated using the LaTeX2HTML translator Version 2K.1beta (1.61)

Copyright © 1993, 1994, 1995, 1996, Nikos Drakos, Computer Based Learning Unit, University of Leeds.
Copyright © 1997, 1998, 1999, Ross Moore, Mathematics Department, Macquarie University, Sydney.

The command line arguments were:
latex2html -split +0 lecture.tex

The translation was initiated by amth247 on 2003-04-07

Footnotes

... errors¹

Heath, Example 2.15, gives a theoretical example to show the importance of pivoting. We will see a practical example in Practical 6.