Next: Further Details: Error Bounds Up: Accuracy and Stability Previous: Further Details: Error Bounds   Contents   Index

# Error Bounds for Linear Least Squares Problems

The linear least squares problem is to find x that minimizes . 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 be the solution computed by one of the driver routines xGELS, xGELSX, xGELSY, xGELSS, or xGELSD (see section 2.3.2). An approximate error bound4.10

may be computed in one of the following ways, depending on which type of driver routine is used:

1.
Suppose the simple driver SGELS is used:

      EPSMCH = SLAMCH( 'E' )
*     Get the 2-norm of the right hand side B
BNORM = SNRM2( M, B, 1 )
*     Solve the least squares problem; the solution X overwrites B
CALL SGELS( 'N', M, N, 1, A, LDA, B, LDB, WORK, LWORK, INFO )
IF ( MIN(M,N) .GT. 0 ) THEN
*        Get the 2-norm of the residual A*X-B
RNORM = SNRM2( M-N, B( N+1 ), 1 )
*        Get the reciprocal condition number RCOND of A
CALL STRCON('I', 'U', 'N', N, A, LDA, RCOND, WORK, IWORK, INFO)
RCOND = MAX( RCOND, EPSMCH )
IF ( BNORM .GT. 0.0 ) THEN
SINT = RNORM / BNORM
ELSE
SINT = 0.0
ENDIF
COST = MAX( SQRT( (1.0E0 - SINT)*(1.0E0 + SINT) ), EPSMCH )
TANT = SINT / COST
ERRBD = EPSMCH*( 2.0E0/(RCOND*COST) + TANT / RCOND**2 )
ENDIF


For example4.11, if ,

then, to 4 decimal places,

, , , , and the true error is .

2.
Suppose the expert driver SGELSX or SGELSY is used. This routine has an input argument RCND, which is used to determine the rank of the input matrix (briefly, the matrix is considered not to have full rank if its condition number exceeds 1/RCND). The code fragment below only computes error bounds if the matrix has been determined to have full rank. When the matrix does not have full rank, computing and interpreting error bounds is more complicated, and the reader is referred to the next section.

      EPSMCH = SLAMCH( 'E' )
*     Get the 2-norm of the right hand side B
BNORM = SNRM2( M, B, 1 )
*     Solve the least squares problem; the solution X overwrites B
RCND = 0
CALL SGELSX( M, N, 1, A, LDA, B, LDB, JPVT, RCND, RANK, WORK,
$INFO ) IF ( RANK.LT.N ) THEN PRINT *,'Matrix less than full rank' ELSE IF ( MIN( M,N ) .GT. 0 ) THEN * Get the 2-norm of the residual A*X-B RNORM = SNRM2( M-N, B( N+1 ), 1 ) * Get the reciprocal condition number RCOND of A CALL STRCON('I', 'U', 'N', N, A, LDA, RCOND, WORK, IWORK, INFO) RCOND = MAX( RCOND, EPSMCH ) IF ( BNORM .GT. 0.0 ) THEN SINT = RNORM / BNORM ELSE SINT = 0.0 ENDIF COST = MAX( SQRT( (1.0E0 - SINT)*(1.0E0 + SINT) ), EPSMCH ) TANT = SINT / COST ERRBD = EPSMCH*( 2.0E0/(RCOND*COST) + TANT / RCOND**2 ) END IF  The numerical results of this code fragment on the above A and b are the same as for the first code fragment. 3. Suppose the other type of expert driver SGELSS or SGELSD is used. This routine also has an input argument RCND, which is used to determine the rank of the matrix A. The same code fragment can be used to compute error bounds as for SGELSX or SGELSY, except that the call to SGELSX must be replaced by:  CALL SGELSD( M, N, 1, A, LDA, B, LDB, S, RCND, RANK, WORK, LWORK,$             IWORK, INFO )


and the call to STRCON must be replaced by:

         RCOND = S( N ) / S( 1 )


Applied to the same A and b as above, the computed is nearly the same, , , and the true error is .

Next: Further Details: Error Bounds Up: Accuracy and Stability Previous: Further Details: Error Bounds   Contents   Index
Susan Blackford
1999-10-01