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

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
bound^{4.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 example

^{4.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 aboveand*A*are the same as for the first code fragment.*b* - 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. 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:*A*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

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