Iterative Refinement - Residual Correction



    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  [Graphics:Images/IterativeRefinementMod_gr_2.gif], is found in three steps:  

1.  Construct the matrices  [Graphics:Images/IterativeRefinementMod_gr_3.gif], if possible.  
2.  Solve  [Graphics:Images/IterativeRefinementMod_gr_4.gif]  for  [Graphics:Images/IterativeRefinementMod_gr_5.gif]  using forward substitution.
3.  Solve  [Graphics:Images/IterativeRefinementMod_gr_6.gif]  for  [Graphics:Images/IterativeRefinementMod_gr_7.gif]  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  [Graphics:Images/IterativeRefinementMod_gr_8.gif].  is an arrangement  [Graphics:Images/IterativeRefinementMod_gr_9.gif]  of these integers in a definite order.  The standard base vectors  [Graphics:Images/IterativeRefinementMod_gr_10.gif]  for  [Graphics:Images/IterativeRefinementMod_gr_11.gif]  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  [Graphics:Images/IterativeRefinementMod_gr_12.gif].  The elements of  [Graphics:Images/IterativeRefinementMod_gr_13.gif]  have the form  


Theorem.  Suppose that  [Graphics:Images/IterativeRefinementMod_gr_15.gif]  is a permutation matrix.  The product  PA  is a new matrix whose rows consists of the rows of  A  rearranged in the new order  [Graphics:Images/IterativeRefinementMod_gr_16.gif].  


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  [Graphics:Images/IterativeRefinementMod_gr_20.gif], is found in four steps:  

1.  Construct the matrices  [Graphics:Images/IterativeRefinementMod_gr_21.gif].  
2.  Compute the column vector   [Graphics:Images/IterativeRefinementMod_gr_22.gif].
3.  Solve  [Graphics:Images/IterativeRefinementMod_gr_23.gif]  for  [Graphics:Images/IterativeRefinementMod_gr_24.gif]  using forward substitution.
4.  Solve  [Graphics:Images/IterativeRefinementMod_gr_25.gif]  for  [Graphics:Images/IterativeRefinementMod_gr_26.gif]  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  [Graphics:Images/IterativeRefinementMod_gr_29.gif].  
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  [Graphics:Images/IterativeRefinementMod_gr_55.gif]  is computed only once, and we have

Use the LU solver to construct  [Graphics:Images/IterativeRefinementMod_gr_57.gif]  which is an approximate solution to  [Graphics:Images/IterativeRefinementMod_gr_58.gif],  i.e.  [Graphics:Images/IterativeRefinementMod_gr_59.gif]   

Form the error in the computation  [Graphics:Images/IterativeRefinementMod_gr_60.gif],  this is called the residual


Now subtract  
[Graphics:Images/IterativeRefinementMod_gr_62.gif]  from  [Graphics:Images/IterativeRefinementMod_gr_63.gif],  obtaining  [Graphics:Images/IterativeRefinementMod_gr_64.gif]  or  [Graphics:Images/IterativeRefinementMod_gr_65.gif].  Use the substitution  [Graphics:Images/IterativeRefinementMod_gr_66.gif]  and this equation becomes  


Use the LU solver to compute  [Graphics:Images/IterativeRefinementMod_gr_68.gif]  as follows  


Then the improvement  [Graphics:Images/IterativeRefinementMod_gr_70.gif]  is made  

(4)            [Graphics:Images/IterativeRefinementMod_gr_71.gif]

The process can be iterated to obtain a sequence [Graphics:Images/IterativeRefinementMod_gr_72.gif]  converging to  [Graphics:Images/IterativeRefinementMod_gr_73.gif]  
    For  [Graphics:Images/IterativeRefinementMod_gr_74.gif]  

Proof  Iterative Refinement


Algorithm (Iterative Refinement).  Compute the factorization  [Graphics:Images/IterativeRefinementMod_gr_76.gif]  is computed once.  The notation  [Graphics:Images/IterativeRefinementMod_gr_77.gif]  means that [Graphics:Images/IterativeRefinementMod_gr_78.gif] is the computed solution for the equation  [Graphics:Images/IterativeRefinementMod_gr_79.gif].  

        Set  [Graphics:Images/IterativeRefinementMod_gr_80.gif]
    For  [Graphics:Images/IterativeRefinementMod_gr_82.gif]  

Computer Programs  Iterative Refinement

Mathematica Subroutine (Iterative Refinement).  The following subroutine will perform the [Graphics:Images/IterativeRefinementMod_gr_84.gif] step of iterative refinement.  Input [Graphics:Images/IterativeRefinementMod_gr_85.gif] and output [Graphics:Images/IterativeRefinementMod_gr_86.gif].  The local variables in the subroutine are labeled as the first step.  The subroutine call should be  [Graphics:Images/IterativeRefinementMod_gr_87.gif],  it is assumed that [Graphics:Images/IterativeRefinementMod_gr_88.gif] is already a global variable and will create and return  [Graphics:Images/IterativeRefinementMod_gr_89.gif].  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  [Graphics:Images/IterativeRefinementMod_gr_92.gif]  of a linear system  [Graphics:Images/IterativeRefinementMod_gr_93.gif].  We see an inherent difficulty in this process:  if we cannot solve  [Graphics:Images/IterativeRefinementMod_gr_94.gif]  exactly, we cannot expect to solve  [Graphics:Images/IterativeRefinementMod_gr_95.gif]  exactly!  The signal feature of iterative improvement is to calculate the residuals  [Graphics:Images/IterativeRefinementMod_gr_96.gif] 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  [Graphics:Images/IterativeRefinementMod_gr_97.gif].  Use iterative refinement.

Solution 2.


Example 3. Solve  [Graphics:Images/IterativeRefinementMod_gr_108.gif].  Use iterative refinement.

Solution 3.


    A linear system with a Hilbert matrix  [Graphics:Images/IterativeRefinementMod_gr_119.gif]  is difficult to solve numerically. The following examples illustrate this situation.  


Example 4.  Solve  [Graphics:Images/IterativeRefinementMod_gr_120.gif],  where [Graphics:Images/IterativeRefinementMod_gr_121.gif] is the Hilbert matrix and  [Graphics:Images/IterativeRefinementMod_gr_122.gif].  Use iterative refinement.

Solution 4.


Example 5.  Solve  [Graphics:Images/IterativeRefinementMod_gr_160.gif],  where [Graphics:Images/IterativeRefinementMod_gr_161.gif] is the Hilbert matrix and  [Graphics:Images/IterativeRefinementMod_gr_162.gif].  Use iterative refinement.
Solution 5.


Example 6.  Solve  [Graphics:Images/IterativeRefinementMod_gr_205.gif],  where [Graphics:Images/IterativeRefinementMod_gr_206.gif] is the Hilbert matrix and  [Graphics:Images/IterativeRefinementMod_gr_207.gif].  Use iterative refinement.
Solution 6.


Example 7.  Solve  [Graphics:Images/IterativeRefinementMod_gr_251.gif],  where [Graphics:Images/IterativeRefinementMod_gr_252.gif] is the Hilbert matrix and  [Graphics:Images/IterativeRefinementMod_gr_253.gif].   Use iterative refinement.
Solution 7.


Example 8.  Solve  [Graphics:Images/IterativeRefinementMod_gr_297.gif],  where [Graphics:Images/IterativeRefinementMod_gr_298.gif] is the Hilbert matrix and  [Graphics:Images/IterativeRefinementMod_gr_299.gif].  Use iterative refinement.
Solution 8.


Example 9.  Solve  [Graphics:Images/IterativeRefinementMod_gr_343.gif],  where [Graphics:Images/IterativeRefinementMod_gr_344.gif] is the Hilbert matrix and  [Graphics:Images/IterativeRefinementMod_gr_345.gif].  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