Error Bounds for Linear Least Squares Problems

The linear least squares problem is to find x that minimizes $\Vert Ax-b \Vert _2$. 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 $\hat{x}$ 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

$$\frac{\Vert \hat{x} - x \Vert _2}{\Vert x \Vert _2} \mathrel... ...box{-.75ex}{$\mathop{\sim}\limits^{\textstyle <}$}}{\tt ERRBD}$$

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 ${\tt SLAMCH('E')} = 2^{-24} = 5.961 \cdot 10^{-8}$,

$$A = \left( \begin{array}{ccc} 4 & 3 & 5 \\ 2 & 5 & 8 \\ 3 & ... ...n{array}{c} 100.1 \\ .1 \\ .01 \\ .01 \end{array} \right) \; ,$$

then, to 4 decimal places,

$$x = \hat{x} = \left( \begin{array}{c} 38.49 \\ 21.59 \\ -23.88 \end{array} \right) \; \; ,$$

${\tt BNORM} = 100.1$, ${\tt RNORM} = 8.843$, ${\tt RCOND} = 4.712 \cdot 10^{-2}$, ${\tt ERRBD} = 4.9 \cdot 10^{-6}$, and the true error is $4.6 \cdot 10^{-7}$.

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 $\hat{x}$ is nearly the same, ${\tt RCOND} = 5.428 \cdot 10^{-2}$, ${\tt ERRBD} = 4.0 \cdot 10^{-6}$, and the true error is $6.6 \cdot 10^{-7}$.