next up previous contents index
Next: Error Bounds for the Up: Accuracy and Stability Previous: Error Bounds for Linear

Error Bounds for Linear Least Squares Problems

 

The linear least squares problem is to find x that minimizes tex2html_wrap_inline18475. We discuss error bounds for the most common case where A is m-by-n with m > n, and A has full rank ; this is called an overdetermined least squares problem   (the following code fragments deal with m=n as well).

Let tex2html_wrap_inline17482 be the solution computed by the driver routine PxGELS (see section 3.2.2). An approximate error bound     
displaymath18463
may be computed in the following way:

      EPSMCH = PSLAMCH( ICTXT, 'E' )
*     Get the 2-norm of the right hand side B
      BNORM = PSLANGE( 'F', M, 1, B, IB, JB, DESCB, WORK )
*     Solve the least squares problem; the solution X overwrites B
      CALL PSGELS( 'N', M, N, 1, A, IA, JA, DESCA, B, IB, JB, DESCB,
     $             WORK, LWORK, INFO )
      IF( MIN( M, N ).GT.0 ) THEN
*        Get the 2-norm of the residual A*X-B
         RNORM = PSLANGE( 'F', M-N, 1, B, IB+N, JB, DESCB, WORK ) 
*        Get the reciprocal condition number RCOND of A
         CALL PSTRCON( 'I', 'U', 'N', N, A, IA, JA, DESCA, RCOND, WORK,
     $                 LWORK, IWORK, LIWORK, INFO )
         RCOND = MAX( RCOND, EPSMCH )
         IF( BNORM.GT.0.0 ) THEN
            SINT = RNORM / BNORM
         ELSE
            SINT = 0.0
         END IF
         COST = MAX( SQRT( ( 1.0E0 - SINT )*( 1.0E0 + SINT ) ), EPSMCH )
         TANT = SINT / COST
         ERRBD = EPSMCH*( 2.0E0 / ( RCOND*COST ) + TANT / RCOND**2 )
      END IF
 

For example, if tex2html_wrap_inline18242,
displaymath18464
then, to four decimal places,
displaymath18465
tex2html_wrap_inline18495, tex2html_wrap_inline18497, tex2html_wrap_inline18499, tex2html_wrap_inline18501, and the true error is tex2html_wrap_inline18503.

Note that in the preceding code fragment, the routine PSLANGE was used to compute the two-norm of the right hand side matrix B and the residual tex2html_wrap_inline18507. This routine was chosen because the result of the computation (BNORM or RNORM, respectively) is automatically known on all process columns within the process grid. The routine PSNRM2 could have also been used to perform this calculation; however, the use of PSNRM2 in this example would have required an additional communication broadcast, because the resulting value of BNORM or RNORM, respectively, is known only within the process column owning B(:,JB).

Further Details: Error Bounds for Linear Least Squares Problems 

The conventional error analysis of linear least squares problems goes as follows . As above, let tex2html_wrap_inline17482 be the solution to minimizing tex2html_wrap_inline18475 computed by ScaLAPACK using the least squares driver PxGELS (see section 3.2.2). We discuss the most common case, where A is overdetermined  (i.e., has more rows than columns) and has full rank [71]:     

The computed solution tex2html_wrap_inline17482 has a small normwise backward error. In other words, tex2html_wrap_inline17482 minimizes tex2html_wrap_inline18519, where E and f satisfy    
displaymath18466
and p(n) is a modestly growing function of n. We take p(n)=1 in the code fragments above. Let tex2html_wrap_inline18531 (approximated by 1/RCOND in the above code fragments), tex2html_wrap_inline18533 (= RNORM above), and tex2html_wrap_inline18535 (SINT = RNORM / BNORM above). Here, tex2html_wrap_inline17638 is the acute angle between the vectors tex2html_wrap_inline18539 and b.     Then when tex2html_wrap_inline18543 is small, the error tex2html_wrap_inline18545 is bounded by
displaymath18467
where tex2html_wrap_inline18549 = COST and tex2html_wrap_inline18551 = TANT in the code fragments above.

We avoid overflow by making sure RCOND and COST are both at least tex2html_wrap_inline18308 EPSMCH, and by handling the case of a zero B matrix separately (BNORM = 0).    

tex2html_wrap_inline18531 may be computed directly from the singular values of A returned by PxGESVD. It may also be approximated by using PxTRCON following calls to PxGELS. PxTRCON estimates tex2html_wrap_inline18559 or tex2html_wrap_inline18561 instead of tex2html_wrap_inline18563, but these can differ from tex2html_wrap_inline18563 by at most a factor of n.          


next up previous contents index
Next: Error Bounds for the Up: Accuracy and Stability Previous: Error Bounds for Linear

Susan Blackford
Tue May 13 09:21:01 EDT 1997