Error Bounds for the Generalized Singular Value Decomposition

The generalized (or quotient) singular value decomposition of an m-by-n matrix A and a p-by-n matrix B is the pair of factorizations

$$A = U \Sigma_1 [0,R] Q^T \; \; {\rm and} \; \; B = V \Sigma_2 [0,R] Q^T$$

where U, V, Q, R, $\Sigma_1$ and $\Sigma_2$ are defined as follows.

• U is m-by-m, V is p-by-p, Q is n-by-n, and all three matrices are orthogonal. If A and B are complex, these matrices are unitary instead of orthogonal, and Q^T should be replaced by Q^H in the pair of factorizations.
• R is r-by-r, upper triangular and nonsingular. [0,R] is r-by-n. The integer r is the rank of $\left( \begin{array}{c} A \\ B \end{array} \right)$, and satisfies $r \leq n$.
• $\Sigma_1$ is m-by-r, $\Sigma_2$ is p-by-r, both are real, nonnegative and diagonal, and $\Sigma_1^T \Sigma_1 + \Sigma_2^T \Sigma_2 = I$. Write $\Sigma_1^T \Sigma_1 = {\rm diag} ( \alpha_1^2 , \ldots , \alpha_r^2 )$ and $\Sigma_2^T \Sigma_2 = {\rm diag} ( \beta_1^2 , \ldots , \beta_r^2 )$, where $\alpha_i$ and $\beta_i$ lie in the interval from 0 to 1. The ratios $\alpha_1 / \beta_1 , \ldots,, \alpha_r / \beta_r$ are called the generalized singular values of the pair A, B. If $\beta_i = 0$, then the generalized singular value $\alpha_i / \beta_i$ is infinite. For details on the structure of $\Sigma_1$, $\Sigma_2$ and R, see section 2.3.5.3.

The generalized singular value decomposition is computed by driver routine xGGSVD (see section 2.3.5.3). We will give error bounds for the generalized singular values in the common case where $\left( \begin{array}{c} A \\ B \end{array} \right)$ has full rank r=n. Let $\hat{\alpha}_i$ and $\hat{\beta}_i$ be the values of $\alpha_i$ and $\beta_i$, respectively, computed by xGGSVD. The approximate error bound^4.10for these values is

$$\vert \hat{\alpha}_i - \alpha_i \vert + \vert \hat{\beta}_i - \beta_i \vert \leq {\tt SERRBD} \; \; .$$

Note that if $\beta_i$ is close to zero, then a true generalized singular value $\alpha_i / \beta_i$ can differ greatly in magnitude from the computed generalized singular value $\hat{\alpha}_i / \hat{\beta}_i$, even if SERRBD is close to its minimum $\epsilon$.

Here is another way to interpret SERRBD: if we think of $\alpha_i$ and $\beta_i$ as representing the subspace $\cal S$ consisting of the straight line through the origin with slope $\alpha_i / \beta_i$, and similarly $\hat{\alpha}_i$ and $\hat{\beta}_i$ representing the subspace $\hat{\cal S}$, then ${\tt SERRBD}$ bounds the acute angle between $\cal S$ and $\hat{\cal S}$. Note that any two lines through the origin with nearly vertical slopes (very large $\alpha / \beta$) are close together in angle. (This is related to the chordal distance in section 4.10.1.)

SERRBD can be computed by the following code fragment, which for simplicity assumes $m \geq n$. (The assumption r=n implies only that $p+m \geq n$. Error bounds can also be computed when $p+m \geq n > m$, with slightly more complicated code.)

      EPSMCH = SLAMCH( 'E' )
*     Compute generalized singular values of A and B
      CALL SGGSVD( 'N', 'N', 'N', M, N, P, K, L, A, LDA, B,
     $             LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ,
     $             WORK, IWORK, INFO )
*     Compute rank of [A',B']'
      RANK = K+L
      IF( INFO.GT.0 ) THEN
         PRINT *,'SGGSVD did not converge'
      ELSE IF( RANK.LT.N ) THEN
         PRINT *,'[A**T,B**T]**T not full rank'
      ELSE IF ( M .GE. N .AND. N .GT. 0 ) THEN
*        Compute reciprocal condition number RCOND of R
         CALL STRCON( 'I', 'U', 'N', N, A, LDA, RCOND, WORK, IWORK,
     $                INFO )
         RCOND = MAX( RCOND, EPSMCH )
         SERRBD = EPSMCH / RCOND
      END IF

For example^4.11, if ${\tt SLAMCH('E')} = 2^{-24} = 5.961 \cdot 10^{-8}$,

$$A = \left( \begin{array}{ccc} 1 & 2 & 3 \\ 3 & 2 & 1 \\ 4 & ... ...egin{array}{ccc} -2 & -3 & 3 \\ 4 & 6 & 5 \end{array} \right)$$

then, to 4 decimal places,

$$\left( \begin{array}{c} \alpha_1 \\ \alpha_2 \\ \alpha_3 \en... ...ft( \begin{array}{c} 0 \\ .6053 \\ .9968 \end{array} \right) ,$$

${\tt SERRBD} = 1.4 \cdot 10^{-6}$, and the true errors are 0, $4.3 \cdot 10^{-7}$ and $1.5 \cdot 10^{-7}$.