\documentclass[11pt]{article}

\usepackage{enumerate}
\usepackage{ifthen}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{color}

%\pagestyle{empty}
\topmargin -0.75in
\oddsidemargin -0.25in
\evensidemargin -0.25in
\textheight 9.25in
\textwidth 7.0in 

\definecolor{ignoreGray}{rgb}{0.75,0.75,0.75}
\newcommand{\ig}{\color{ignoreGray}}

\newcommand{\myspace}{\mbox{\parbox[c][0.35in]{0.35in}{\hfill}}}
\newcommand{\varspace}[2]{\mbox{\parbox[c][#1]{#2}{\hfill}}}

\newcommand{\B}[1]{{\bf #1}}

\newcommand{\E}[2]{\B e_{#1}\B e_{#2}^T}

\newcommand{\force}{{$\mbox{\;}$}}

\begin{document}

\noindent \hfill
\fbox{\fbox{\parbox{6.9in}{{\scshape{Linear Algebra \hfill Partial Pivoting \& PLU-Decomposition}} }}} \hfill 

\vspace{0.2in}

Gauss-Jordan elimination is the work horse of linear algebra. While it gets the job done in theory, it is not stable {\it numerically}. If you are working with approximate coefficients in your linear system, regular Gauss-Jordan elimination can significantly amplify error. To combat this we have the methods of partial pivoting and full pivoting. 

Partial pivoting is a minor modification of Gauss-Jordan elimination. With regular elimination, after locating a pivot position, you merely swap rows so there is something non-zero in the pivot position. Partial pivoting demands that you always swap the largest (in magnitude) possible entry into the pivot position. 

Why does this help? Well, in the forward pass, error is created when we add a {\it multiple} of a row to another row. The multiple we use can magnify whatever round off error is already there. By choosing the largest available pivot, we guarantee that our multiples always have magnitudes of 1 or less. This means we limit the propagation of error (possibly significantly). 

However, this does not always fix the instability issue. In theory, there are systems that (when solved with partial pivoting) small perturbations of inputs yields large perturbations of outputs. On the other hand, the technique of {\it full pivoting} is stable. Full pivoting means you swap rows {\it and columns} if necessary to get a pivot of largest magnitude. While full pivoting is not hard to implement, these column swaps are quite annoying. This amounts to renaming variables. So when you finish row reduction, you then have to track down how your variables where permuted. In many circumstances, partial pivoting is good enough. 

\vspace{0.15in}

\noindent {\bf\large Example:} $\begin{array}{rcccccl} x & + & 4y & - & 2z & = & 3 \\ -x & + & 3y & + & z & = & 8 \\ 2x & & & + & 3z & = & 11 \end{array}$. We have $A\B x = \B b$ where $A=\begin{bmatrix} 1 & 4 & -2 \\ -1 & 3 & 1 \\ 2 & 0 & 3 \end{bmatrix}$ and $\B b=\begin{bmatrix} 3 \\ 8 \\ 11 \end{bmatrix}$.

We will solve this system by doing a partial pivoting forward pass and back substitution. 

\vspace{-0.1in}

$$\hspace*{-0.25in} [A\; : \; \B b] = \begin{bmatrix} 1 & 4 & -2 & : & 3 \\ -1 & 3 & 1 & : & 8 \\ 2 & 0 & 3 & : & 11\end{bmatrix} 
\quad \stackrel{R_1 \leftrightarrow R_3}{\widetilde{\qquad}} \quad
\begin{bmatrix} 2 & 0 & 3 & : & 11 \\ -1 & 3 & 1 & : & 8 \\ 1 & 4 & -2 & : & 3 \end{bmatrix} 
\quad \stackrel{\frac12 R_1 + R_2}{\widetilde{\qquad}}
\quad \stackrel{-\frac12 R_1 + R_3}{\widetilde{\qquad}} \quad
\begin{bmatrix} 2 & 0 & 3 & : & 11 \\ 0 & 3 & \frac52  & : & \frac{27}2 \Large\strut \\ 0 & 4 & -\frac72 & : & -\frac52 \Large\strut \end{bmatrix} 
$$

$$\quad \stackrel{R_2 \leftrightarrow R_3}{\widetilde{\qquad}} \quad
\begin{bmatrix} 2 & 0 & 3 & : & 11 \\ 0 & 4 & -\frac72 & : & -\frac52 \Large\strut \\ 0 & 3 & \frac52  & : & \frac{27}2 \Large\strut \end{bmatrix} 
\quad \stackrel{-\frac34R_2 + R_3}{\widetilde{\qquad}}
\begin{bmatrix} 2 & 0 & 3 & : & 11 \\ 0 & 4 & -\frac72 & : & -\frac52 \Large\strut \\ 0 & 0 & \frac{41}8  & : & \frac{123}8 \Large\strut \end{bmatrix} 
$$

Notice our row swaps to make sure we had pivots of maximal magnitudes. This completes our partial pivoting forward pass. Now we back substitute to solve. 

The last row of our augmented matrix says $\frac{41}8z=\frac{123}8$ and so $z=3$. The second row says $4y-\frac72z=-\frac52$ and so $y=\frac14\left(-\frac52+\frac72(3)\right)=2$. Finally, the top row says $2x+3y=11$ and so $x=\frac12\left(11-3(3)\right)=1$. 

Therefore, the (unique) solution is $x=1$, $y=2$, and $z=3$. Of course, it would have been much easier to just use standard row reduction!

\vspace{0.15in}

Next, we will compute a PLU-decomposition of $A$ using partial pivoting. Then we will use this decomposition to once again solve our system. To compute a PLU-decomposition we need to complete a forward pass on $A$. This time we will keep track of the type III operation multiples that show up. These will be denoted by {\color{red}red entries} in parentheses. In particular, if we perform ``$a R_i + R_j$'', then we will put ``${\color{red}(-a)}$'' in the $j$-th row directly below our pivot (i.e., the entry that was being wiped out by this operation).

Because of the way that row swaps and these type III operations interact, we want to move these entries along with the others when we swap rows. Following this procedure will leave us with entries that go below the diagonal in our lower triangular matrix $L$. In other words, we will be able to read off $L$ and $U$ at the end of the forward pass. For the permutation piece, it turns out our permutations could have been done up front, so we end up with $(\text{Permute})A = LU$. Thus the permutations done actually build $P^{-1}$. We have $P^{-1}A=LU$ so that $A=PLU$. Therefore, the $P$ we're looking for is built by {\it undoing} the swaps done in the forward pass (i.e., take the identity matrix $I$ and then do the swaps from last to first -- this gives us our $P$).

$$A = \begin{bmatrix} 1 & 4 & -2 \\ -1 & 3 & 1 \\ 2 & 0 & 3 \end{bmatrix} 
\quad \stackrel{R_1 \leftrightarrow R_3}{\widetilde{\qquad}} \quad
\begin{bmatrix} 2 & 0 & 3 \\ -1 & 3 & 1 \\ 1 & 4 & -2 \end{bmatrix} 
\quad \stackrel{\frac12 R_1 + R_2}{\widetilde{\qquad}}
\begin{bmatrix} 2 & 0 & 3 \\ {\color{red}(-\frac12)} & 3 & \frac52  \Large\strut \\ 1 & 4 & -2 \end{bmatrix} 
\quad \stackrel{-\frac12 R_1 + R_3}{\widetilde{\qquad}} \quad
\begin{bmatrix} 2 & 0 & 3 \\ {\color{red}(-\frac12)} & 3 & \frac52  \Large\strut \\ {\color{red}(\frac12)} & 4 & -\frac72 \Large\strut \end{bmatrix} 
$$

$$\quad \stackrel{R_2 \leftrightarrow R_3}{\widetilde{\qquad}} \quad
\begin{bmatrix} 2 & 0 & 3 \\ {\color{red}(\frac12)} & 4 & -\frac72 \Large\strut \\ {\color{red}(-\frac12)} & 3 & \frac52  \Large\strut \end{bmatrix} 
\quad \stackrel{-\frac34R_2 + R_3}{\widetilde{\qquad}}
\begin{bmatrix} 2 & 0 & 3 \\ {\color{red}(\frac12)} & 4 & -\frac72 \Large\strut \\ {\color{red}(-\frac12)} & {\color{red}(\frac34)} & \frac{41}8  \Large\strut \end{bmatrix} 
$$

$$I 
\  \stackrel{R_1 \leftrightarrow R_3}{\widetilde{\qquad}} \ 
\begin{bmatrix} 0 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 & 0 \end{bmatrix}
\ \stackrel{R_2 \leftrightarrow R_3}{\widetilde{\qquad}} \ 
\begin{bmatrix} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix} = P^{-1}
\qquad
I
\  \stackrel{R_2 \leftrightarrow R_3}{\widetilde{\qquad}} \ 
\begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \end{bmatrix}
\  \stackrel{R_1 \leftrightarrow R_3}{\widetilde{\qquad}} \ 
\begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{bmatrix}=P
$$

Therefore, $A = \begin{bmatrix} 1 & 4 & -2 \\ -1 & 3 & 1 \\ 2 & 0 & 3 \end{bmatrix} = PLU = \underbrace{\begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{bmatrix}}_P
\underbrace{\begin{bmatrix} 1 & 0 & 0 \\ \frac12 & 1 & 0 \Large\strut \\ -\frac12 & \frac34 & 1 \Large\strut \end{bmatrix}}_L
\underbrace{\begin{bmatrix} 2 & 0 & 3 \\ 0 & 4 & -\frac72 \Large\strut \\ 0 & 0 & \frac{41}8 \Large\strut \end{bmatrix}}_U
$

To use this decomposition to solve $A\B x=\B b$ we first take care of the permutation: $PLU\B x = \B b$ so $LU\B x = P^{-1}\B b$. In other words, we need to permute $\B b$ according to the swaps done in our forward pass (from first to last). Alternatively, we could just multiply by $P^{-1}$ (as computed above).
$$\B b = \begin{bmatrix} 3 \\ 8 \\ 11 \end{bmatrix} 
\quad \stackrel{R_1 \leftrightarrow R_3}{\widetilde{\qquad}} \quad
\begin{bmatrix} 11 \\ 8 \\ 3 \end{bmatrix}
\quad \stackrel{R_2 \leftrightarrow R_3}{\widetilde{\qquad}} \quad
\begin{bmatrix} 11 \\ 3 \\ 8 \end{bmatrix} = \B c$$

Next, let $U\B x=\B y$ so that $LU\B x=P^{-1}\B b$ is $L\B y = \B c$. We solve this system by forward substitution. Notice that the first row of $[L\; :\; \B c]$ would tell us that $y_1 =11$. The next row tells us $\frac12y_1+y_2=3$ so that $y_2=3-\frac12(11)=-\frac52$. The final row says $-\frac12y_1+\frac34y_2+y_3=8$ so that $y_3=8+\frac12(11)-\frac34(-\frac52)=\frac{123}8$ so $\B y = \left[11\;-\frac52\;\;\frac{123}8\right]^T$.

Finally, we solve $U\B x=\B y$ using forward substitution. The last row of $[U\;:\;\B y]$ would tell us that $\frac{41}8x_3=\frac{123}8$ so $x_3=3$. The middle row tells us that 
$4x_2-\frac72x_3=-\frac52$ so that $x_2 =\frac14\left(-\frac52+\frac72(3)\right) = 2$. Last of all, the first row says $2x_1+3x_3=11$ so that $x_1=\frac12(11-3(3))=1$. Therefore, once again, we find that $x=1$, $y=2$, and $z=3$ is the unique solution of our system.

\vspace{0.2in}

Of course all of this is a lot of work to solve a fairly simple system. So why? We already discussed that partial pivoting helps with numerical issues. The PLU-decomposition is helpful if are solving multiple systems with the same coefficient matrix $A$. If so, you find $A=PLU$ once and then use it over and over again. With this decomposition in hand, you simply (1) permute the entries of $\B b$, (2) solve a lower triangular system with forward substitution, and (3) solve an upper triangular system with back substitution. These 3 steps can be done {\it very} quickly (even for extremely large systems) and they won't introduce much new round off error (most of that is generated when computing the PLU-decomposition). 


\end{document}
