Error Bounds for Linear Equation Solving

 

Let Ax=b be the system to be solved, and the computed solution. Let n be the dimension of A. An approximate error bound  for may be obtained in one of the following two ways, depending on whether the solution is computed by a simple driver or an expert driver:

1. Suppose that Ax=b is solved using the simple driver PSGESV    (section 3.2.1). Then the approximate error bound
   
   
   
   can be computed by the following code fragment.

         EPSMCH = PSLAMCH( ICTXT, 'E' )
   *     Get infinity-norm of A
         ANORM = PSLANGE( 'I', N, N, A, IA, JA, DESCA, WORK )
   *     Solve system; The solution X overwrites B
         CALL PSGESV( N, 1, A, IA, JA, DESCA, IPIV, B, IB, JB, DESCB, INFO )
         IF( INFO.GT.0 ) THEN
            PRINT *,'Singular Matrix'
         ELSE IF( N.GT.0 ) THEN
   *        Get reciprocal condition number RCOND of A
            CALL PSGECON( 'I', N, A, IA, JA, DESCA, ANORM, RCOND, WORK,
        $                 LWORK, IWORK, LIWORK, INFO )
            RCOND = MAX( RCOND, EPSMCH )
            ERRBD = EPSMCH / RCOND
         END IF

    

   For example, suppose
   
   ,
   
   Then (to 4 decimal places)
   
   , , the true reciprocal condition number , , and the true error .  

2. Suppose that Ax=b is solved using the expert driver PSGESVX (section 3.2.1).    This routine provides an explicit error bound FERR, measured with the infinity-norm:  
   
   For example, the following code fragment solves Ax=b and computes an approximate error bound FERR:

         CALL PSGESVX( 'E', 'N', N, 1, A, IA, JA, DESCA, AF, IAF, JAF,
        $              DESCAF, IPIV, EQUED, R, C, B, IB, JB, DESCB, X, IX,
        $              JX, DESCX, RCOND, FERR, BERR, WORK, LWORK, IWORK,
        $              LIWORK, INFO )
         IF( INFO.GT.0 ) PRINT *,'(Nearly) Singular Matrix'

   For the same A and b as above, , , and the actual error is .

This example illustrates that the expert driver provides an error bound with less programming effort than the simple driver, and also that it may produce a significantly more accurate answer.

Similar code fragments, with obvious adaptations, may be used with all the driver routines PxPOSV and PxPOSVX      in Table 3.2. For example, if a symmetric positive definite or Hermitian positive definite system is solved by using the simple driver PxPOSV,      then PxLANSY or PxLANHE, respectively, must be used to compute ANORM, and PxPOCON     must be used to compute RCOND.

The drivers PxGBSV      (for solving general band matrices with partial pivoting), PxPBSV      (for solving positive definite band matrices) and PxPTSV      (for solving positive definite tridiagonal matrices), do not yet have the corresponding routines needed to compute error bounds, namely, PxLAnHE to compute ANORM and PxyyCON to compute RCOND.

The drivers PxDBSV      (for solving general band matrices) and PxDTSV      (for solving general tridiagonal matrices) do not pivot for numerical stability, and so may be faster but less accurate than their pivoting counterparts above. These routines may be used safely when any diagonal pivot sequence leads to a stable factorization; diagonally dominant matrices and symmetric positive definite matrices [71] have this property, for example.

Further Details: Error Bounds for Linear Equation Solving 

The conventional error analysis of linear equation  solving goes as follows. Let Ax=b be the system to be solved. Let be the solution computed by ScaLAPACK (or LAPACK) using any of their linear equation solvers. Let r be the residual . In the absence of rounding error, r would be zero and would equal x; with rounding error, one can only say the following:

The normwise backward error of the computed solution ,     with respect to the infinity norm, is the pair E,f, which minimizes

subject to the constraint . The minimal value of is given by

One can show that the computed solution satisfies , where p(n) is a modestly growing function of n. The corresponding condition number is .   The error is bounded by

In the first code fragment in the preceding section, , which is in the numerical example, is approximated by . Approximations  of -- or, strictly speaking, its reciprocal RCOND -- are returned by computational routines PxyyCON (section 3.3.1) or driver routines                PxyySVX (section 3.2.1). The code fragment makes sure RCOND is at least EPSMCH to avoid overflow in computing ERRBD.   This limits ERRBD to a maximum of 1, which is no loss of generality because a relative error of 1 or more indicates the same thing:    a complete loss of accuracy.   Note that the value of RCOND returned by PxyySVX may apply to a linear system obtained from Ax=b by equilibration, namely, scaling the rows and columns of A in order to make the condition number smaller. This is the case in the second code fragment in the preceding section, where the program chose to scale the rows by the factors returned in , resulting in .

As stated in section 6.4.2, this approach does not respect the presence of zero or tiny entries in A. In contrast, the ScaLAPACK computational routines PxyyRFS (section 3.3.1) or driver routines PxyySVX (section 3.2.1) will (except in rare cases) compute a solution with the following properties:

The componentwise backward error of the computed solution is the pair E,f which minimizes    

(where we interpret 0/0 as 0) subject to the constraint . The minimal value of is given by

One can show that for most problems the computed by PxyySVX satisfies , where p(n) is a modestly growing function of n. In other words, is the exact solution of the perturbed problem , where E and f are small relative perturbations in each entry of A and b, respectively. The corresponding condition number is .   The error is bounded by

The routines PxyyRFS and PxyySVX return ,                              which is called BERR  (for Backward ERRor), and a bound on the the actual error , called FERR   (for Forward ERRor), as in the second code fragment in the last section. FERR is actually calculated by the following formula, which can be smaller than the bound given above:

Here, is the computed value of the residual , and the norm in the numerator is estimated by using the same estimation subroutine used for RCOND.

The value of BERR for the example in the preceding section is .

Even in the rare cases where PxyyRFS fails to make BERR close to its minimum , the error bound FERR may remain small. See [9] for details.