Gauss elimination vs LU factorization

Question

I wanted to know which of the two methods gives more accurate results when solving linear equations problems, Gauss elimination or LU factorization? And what happens if we use partial pivoting while applying the first method, is this more accurate than the LU factorization?

Answer 1

I think the book Numerical Linear Algebra by Trefethen and Bau has a pretty good explanation of this topic. I've summarized some of the content from chapters 20 and 21.

Gaussian Elimination is unstable

Consider the matrix $$ A = \begin{bmatrix} 10^{-20} & 1 \\ 1 & 1\end{bmatrix} $$

This has exact LU decomposition, $$ L = \begin{bmatrix} 1 & 0 \\ 10^{20} & 1\end{bmatrix},~~~~ U = \begin{bmatrix} 10^{-20} & 1 \\ 0 & 1-10^{20}\end{bmatrix} $$

However, suppose we make small (relative) rounding errors and end up representing $1-10^{20}$ as $-10^{20}$. Then, $$ \begin{bmatrix} 1 & 0 \\ 10^{20} & 1\end{bmatrix}\begin{bmatrix} 10^{-20} & 1 \\ 0 & -10^{20}\end{bmatrix} = \begin{bmatrix} 10^{-20} & 1 \\ 1 & 0 \end{bmatrix} $$ which is far from $A$.

LU is Gaussian Elimination without Pivoting

To expand on my comment, when you do Gaussian Elimination without pivoting, each step can be written in terms of matrix multiplication on the left by a lower triangular matrix. The end result is an upper triangular matrix so written in matrix form Gaussian elimination looks something like $$ L_{m-1}\cdots L_2L_1 A = U $$

The product of lower triangular matrices is still lower triangular, as is the inverse, so setting $L^{-1} = L_{m-1}\cdots L_2L_1$ gives the LU decomposition.

As a concrete example, let's consider the matrix $$ A = \begin{bmatrix} 2 & 1 & 1\\ 4 & 3 & 2\\ 8 & 7 & 9 \end{bmatrix} $$

The first step of Gaussian elimination would be to subtract twice the first row from the second, and 4 times the first from the third. We can write this as, $$ L_1 = \begin{bmatrix} 1 \\ -2 & 1\\ -4 & & 1 \end{bmatrix} $$ You can verify that, $$ L_1A = \begin{bmatrix} 2 & 1 & 1 \\ 0 & 1 & 0 \\ 0 & 3 & 5 \end{bmatrix} $$

The next step of Gaussian elimination is to subtract 3 times the second row from the third row. Again, we can write this as, $$ L_2 = \begin{bmatrix} 1 \\ & 1 \\ & -3 & 1 \end{bmatrix} $$ and verify that, $$ L_2(L_1A) = \begin{bmatrix} 2 & 1 & 1 \\ 0 & 1 & 0 \\ 0 & 0 & 5 \end{bmatrix} $$

It turns out that the invese of $L_i$ is quite simple, you just negate the entries below the main diagonal. Similarly, the product of such matrices is also simple, you just add the entries below the diagonal. You can again verify that, $$ L = (L_2L_1)^{-1} = L_1^{-1}L_2^{-1} = \begin{bmatrix} 1 \\ 2 & 1\\ 4 & & 1 \end{bmatrix} \begin{bmatrix} 1 \\ & 1 \\ & 3 & 1 \end{bmatrix} = \begin{bmatrix} 1 \\ 2 & 1 \\ 4 & 3 & 1 \end{bmatrix} $$

PLU is Gaussian Elimination with Pivoting

Now, if you pivot at each step, you can view this as swapping the rows of the working matrix. This is the same as left multiplication by a permutation matrix. So doing gaussian elimination with pivoting will give a set of operations like $$ L_{m-1}P_{m-1}\cdots L_2P_2L_1P_1 A = U $$

If we define, $$ L_k' = P_{m-1}\cdots P_{k+1}L_kP_{k+1}^{-1}\cdots P_{m-1}^{-1} $$ you will see we can write the above series of operations as $$ (L_{m-1}'\cdots L_2'L_1')(P_{m-1}^{-1}\cdots P_1^{-1}) A = U $$

Which you can rearrange to get a PLU decomposition.