Background

    In the LU-Factorization module we developed the "LU-solve" method which will now be used as the computational engine in the iterative refinement - residual correction method.  

Definition (LU-Factorization).  The matrix A has an LU-factorization if it can be expressed as the product of a lower-triangular matrix L and an upper triangular matrix U  

        .

Theorem (A = LU;  Factorization with NO Pivoting).  Assume that A has a LU-factorization.  The solution X to the linear system  , is found in three steps:  

    1.  Construct the matrices  , if possible.  
    2.  Solve    for    using forward substitution.
    3.  Solve    for    using back substitution.

Proof  LU Factorization  LU Factorization  

    The above theorem assumes that there are no row interchanges.  If row interchanges are permitted then a factorization of a "permuted matrix" will be obtained.  A permutation of the first n positive integers  .  is an arrangement    of these integers in a definite order.  The standard base vectors    for    are used in the next definition.  

Definition (Permutation Matrix).  An n×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 and  P  can be written as  .  The elements of    have the form  

          

Theorem.  Suppose that    is a permutation matrix.  The product  PA  is a new matrix whose rows consists of the rows of  A  rearranged in the new order  .  

Exploration

Theorem.  If  A  is a nonsingular matrix, then there exists a permutation matrix  P  so that  PA  has an LU-factorization

        P A  =  L U.  

Theorem (PA = LU;  Factorization with Pivoting).  Given that A is nonsingular.  The solution X to the linear system  , is found in four steps:  

    1.  Construct the matrices  .  
    2.  Compute the column vector   .
    3.  Solve    for    using forward substitution.
    4.  Solve    for    using back substitution.

Proof  LU Factorization  LU Factorization  

Computer Programs  LU Factorization  LU Factorization  

The "LU-solve" Method

    The following pair of subroutines are the heart of the computational engine for the computational engine in the iterative refinement - residual correction method.  When we refer to a "LU-solve step" it means use LUfactor followed by SolveLU.  The next example will review these concepts, and more explanation is in the LU-Factorization module

Mathematica Subroutine (LUfactor).  Matrix  A  is a global variable and elements are changed when the LUfactor is executed.  We save a copy of A in A0.

Mathematica Subroutine (SolveLU).  The subroutine SolveLU which is similar to the back substitution subroutine.

Example 1.  Use the P A  =  L U factorization with pivoting to solve the linear system  .  
Solution 1.

Iterative Refinement - Residual Correction Method

    This is a method for improving the accuracy of a solution produced using LU-factorization.  To start the process, the factorization    is computed only once, and we have
    
(1)            .  

Use the LU solver to construct    which is an approximate solution to  ,  i.e.     

Form the error in the computation  ,  this is called the residual

(2)            .  

Now subtract    from  ,  obtaining    or  .  Use the substitution    and this equation becomes  

               .   

Use the LU solver to compute    as follows  

(3)            .  

Then the improvement    is made  

(4)            

The process can be iterated to obtain a sequence   converging to    
    
            For    
    
            

Proof  Iterative Refinement

Algorithm (Iterative Refinement).  Compute the factorization    is computed once.  The notation    means that is the computed solution for the equation  .  

        Set  
    
        
    
        For    
    
        

Computer Programs  Iterative Refinement

Mathematica Subroutine (Iterative Refinement).  The following subroutine will perform the step of iterative refinement.  Input and output .  The local variables in the subroutine are labeled as the first step.  The subroutine call should be  ,  it is assumed that is already a global variable and will create and return  .  Notice that the computational part of the subroutine consists of only three lines:

        

The remainder of the subroutine consists of print statements.

Remark 1.  Historically, iterative improvement was devised to use a limited amount of "double precision" arithmetic to try to refine the computed solution    of a linear system  .  We see an inherent difficulty in this process:  if we cannot solve    exactly, we cannot expect to solve    exactly!  The signal feature of iterative improvement is to calculate the residuals   in double precision and everything else in single precision.  FORTRAN is well suited for defining the type of a variable to be either single precision or double precision.  A resource for this technique is found in, Section 2.5 Iterative Improvement of a Solution to Linear Equations pp. 47-50, in the book by William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery: Numerical Recipes in Fortran 77, Cambridge University Press, 1992, Cambridge, UK.

Remark 2.  The examples we present are for pedagogical purposes and because this module does not use FORTRAN, we cannot illustrate the full purpose of iterative improvement.  Also, because the convergence will be rapid, we only perform two iterations and to not need a subroutine for the computations.

Example 2.  Solve  .  Use iterative refinement.

Solution 2.

Example 3. Solve  .  Use iterative refinement.

Solution 3.

    A linear system with a Hilbert matrix    is difficult to solve numerically. The following examples illustrate this situation.  

Example 4.  Solve  ,  where is the Hilbert matrix and  .  Use iterative refinement.

Solution 4.

Example 5.  Solve  ,  where is the Hilbert matrix and  .  Use iterative refinement.
Solution 5.

Example 6.  Solve  ,  where is the Hilbert matrix and  .  Use iterative refinement.
Solution 6.

Example 7.  Solve  ,  where is the Hilbert matrix and  .   Use iterative refinement.
Solution 7.

Example 8.  Solve  ,  where is the Hilbert matrix and  .  Use iterative refinement.
Solution 8.

Example 9.  Solve  ,  where is the Hilbert matrix and  .  Use iterative refinement.
Solution 9.

Research Experience for Undergraduates

Iterative Refinement  Internet hyperlinks to web sites and a bibliography of articles.  

Hilbert Matrix  Internet hyperlinks to web sites and a bibliography of articles.  

Ill-Conditioned Linear Systems  Internet hyperlinks to web sites and a bibliography of articles.  

Download this Mathematica Notebook Iterative Refinement

(c) John H. Mathews 2005