Lu Decomposition With Partial Pivoting Example

LU Decomposition with Partial Pivoting: A Comprehensive Guide with Examples

LU decomposition is a fundamental technique in linear algebra used to solve systems of linear equations, find the determinant of a matrix, and compute the inverse of a matrix. It involves factoring a square matrix into a lower triangular matrix (L) and an upper triangular matrix (U). While straightforward for many matrices, incorporating partial pivoting significantly enhances the numerical stability of the algorithm, especially when dealing with ill-conditioned matrices. This article provides a comprehensive explanation of LU decomposition with partial pivoting, including detailed examples and practical considerations.

Understanding LU Decomposition

The core idea behind LU decomposition is to represent a square matrix A as the product of two triangular matrices:

A = LU

where:

• L is a lower triangular matrix with ones on its main diagonal.
• U is an upper triangular matrix.

This factorization simplifies solving the linear system Ax = b because it transforms the problem into two simpler triangular systems:

1. Ly = b (forward substitution)
2. Ux = y (backward substitution)

Solving these triangular systems is computationally much less expensive than directly solving the original system.

The Need for Partial Pivoting

While LU decomposition is efficient, it can be numerically unstable, particularly when dealing with matrices containing small pivot elements (diagonal elements during the Gaussian elimination process). Small pivots can lead to significant round-off errors, producing inaccurate results. Partial pivoting addresses this instability by strategically reordering the rows of the matrix to ensure that the pivot elements are as large as possible in magnitude.

Partial Pivoting Algorithm

The partial pivoting algorithm modifies the standard Gaussian elimination process by incorporating row swaps. At each step of the elimination, the algorithm identifies the element with the largest absolute value in the current column (below the diagonal) and swaps its row with the pivot row. This ensures that the pivot element is always the largest in its column, minimizing the effect of round-off errors.

Here's a step-by-step outline of the algorithm:

1. Initialization: Start with the given matrix A.
2. Iteration: For each column k from 1 to n-1:
   • Find the pivot: Locate the row index i containing the largest absolute value element in column k, starting from row k.
   • Row Swap: If i ≠ k, swap rows i and k.
   • Elimination: For each row j from k+1 to n:
     • Calculate the multiplier: m_jk = a_jk / a_kk
     • Subtract the multiple of row k from row j: a_j = a_j - m_jk * a_k.
3. Result: The resulting matrix is the upper triangular matrix U. The lower triangular matrix L is constructed from the multipliers m_jk. The row swaps are recorded using a permutation matrix P.

Example: LU Decomposition with Partial Pivoting

Let's illustrate LU decomposition with partial pivoting using a concrete example. Consider the following matrix:

A =  [ 2   1   1 ]
     [ 4   3  -1 ]
     [ 8   5   1 ]

Step 1: Find the pivot in the first column

The largest element in the first column is 8 (in the third row). Therefore, we swap the first and third rows:

A' =  [ 8   5   1 ]
      [ 4   3  -1 ]
      [ 2   1   1 ]

Step 2: Elimination

Now we perform Gaussian elimination.

• Row 2: We subtract (4/8) * Row 1 from Row 2:

  [ 8   5   1 ]
  [ 0   1  -3/2]
  [ 2   1   1 ]
  

• Row 3: We subtract (2/8) * Row 1 from Row 3:

  [ 8   5   1 ]
  [ 0   1  -3/2]
  [ 0  -3/4  3/4 ]
  

Step 3: Find the pivot in the second column

The largest element in the second column (below the diagonal) is 1. No row swap is needed.

Step 4: Elimination

• Row 3: We subtract (-3/4) * Row 2 from Row 3:

  [ 8   5   1 ]
  [ 0   1  -3/2]
  [ 0   0  -3/8 ]
  

The resulting upper triangular matrix U is:

U =  [ 8   5   1 ]
     [ 0   1  -3/2]
     [ 0   0  -3/8 ]

The lower triangular matrix L is constructed from the multipliers:

L =  [ 1   0   0 ]
     [ 1/2 1   0 ]
     [ 1/4 3/4 1 ]

The permutation matrix P represents the row swaps:

P =  [ 0   0   1 ]
     [ 0   1   0 ]
     [ 1   0   0 ]

Therefore, the LU decomposition with partial pivoting is: PA = LU.

Solving Linear Systems using LU Decomposition with Partial Pivoting

Once we have the LU decomposition, solving the linear system Ax = b becomes straightforward:

1. Permute: Pb = b' (apply the row swaps from P to b)
2. Forward Substitution: Ly = b'
3. Backward Substitution: Ux = y

Let's say we want to solve Ax = b, where b = [5, 2, 1]^T

1. Pb = b':

   b' = [1, 2, 5]T
   

2. Forward Substitution (Ly = b'):

   y1 = 1
   (1/2)y1 + y2 = 2  => y2 = 3/2
   (1/4)y1 + (3/4)y2 + y3 = 5 => y3 = 11/8
   

3. Backward Substitution (Ux = y):

   -3/8x3 = 11/8 => x3 = -11/3
   x2 - (3/2)x3 = 3/2 => x2 = -2/3
   8x1 + 5x2 + x3 = 1 => x1 = 16/24 = 2/3
   

Therefore, the solution to Ax = b is x = [2/3, -2/3, -11/3]^T

Advantages of Partial Pivoting

• Improved Numerical Stability: Partial pivoting significantly reduces the risk of round-off errors, leading to more accurate solutions, especially for ill-conditioned matrices.
• Avoids Division by Zero: By selecting the largest pivot, partial pivoting prevents division by zero, a common problem in Gaussian elimination.

Disadvantages of Partial Pivoting

• Increased Computational Cost: Partial pivoting adds overhead to the LU decomposition process due to the row swapping and searching for the maximum pivot.
• Not Always Necessary: For well-conditioned matrices, partial pivoting may not be necessary, and the extra computational cost might be unnecessary.

Conclusion

LU decomposition with partial pivoting is a powerful and robust technique for solving systems of linear equations. While the addition of partial pivoting adds computational overhead, the significant improvement in numerical stability makes it an essential tool for many applications, particularly those involving potentially ill-conditioned matrices. Understanding the algorithm and its nuances is crucial for anyone working with numerical linear algebra. This article has provided a detailed explanation with a worked example to solidify the understanding of this important concept. Remember to choose the right numerical method based on the characteristics of your matrix and the desired accuracy of your results.