Further Details: Error Bounds for the Generalized Singular Value Decomposition

The GSVD algorithm used in LAPACK ([83,10,8]) is backward stable:

Let the computed GSVD of A and B be $\hat{U} \hat{\Sigma}_1 [0,\hat{R}] \hat{Q}^T$ and $\hat{V} \hat{\Sigma}_2 [0,\hat{R}] \hat{Q}^T$. This is nearly the exact GSVD of A+E and B+F in the following sense. E and F are small:

$$\Vert E\Vert _2 / \Vert A\Vert _2 \leq p(n) \epsilon \; \; {... ...; \; \Vert F\Vert _2 / \Vert B\Vert _2 \leq p(n) \epsilon \; ;$$

there exist small $\delta \hat{Q}$, $\delta \hat{U}$, and $\delta \hat{V}$ such that $\hat{Q} + \delta \hat{Q}$, $\hat{U}+ \delta \hat{U}$, and $\hat{V}+ \delta \hat{V}$ are exactly orthogonal (or unitary):

$$\Vert \delta \hat{Q} \Vert _2 \leq p(n) \epsilon \; \; , \; ... ...\; \; \Vert \delta \hat{V} \Vert _2 \leq p(n) \epsilon \; \; ;$$

and

$$A + E = ( \hat{U} + \delta \hat{U} ) \hat{\Sigma}_1 [ 0, \ha... ...V} ) \hat{\Sigma}_2 [ 0, \hat{R}] (\hat{Q} + \delta \hat{Q})^T$$

is the exact GSVD of A+E and B+F. Here p(n) is a modestly growing function of n, and we take p(n)=1 in the above code fragment.

Let $\alpha_i$ and $\beta_i$ be the square roots of the diagonal entries of the exact $\Sigma_1^T \Sigma_1$ and $\Sigma_2^T \Sigma_2$, and let $\hat{\alpha}_i$ and $\hat{\beta}_i$ the square roots of the diagonal entries of the computed $\hat{\Sigma}_1^T \hat{\Sigma}_1$ and $\hat{\Sigma}_2^T \hat{\Sigma}_2$. Let

$$G = \left( \begin{array}{c} A \\ B \end{array} \right) \; \;... ...( \begin{array}{c} \hat{A} \\ \hat{B} \end{array} \right) \; .$$

Then provided G and $\hat{G}$ have full rank n, one can show [96,82] that

$$\left( \sum_{i=1}^n [( \hat{\alpha}_i - \alpha_i )^2 + ( \h... ...{\min ( \sigma_{\min} (G) , \sigma_{\min} (\hat{G}) )} \; \; .$$

In the code fragment we approximate the numerator of the last expression by $\epsilon \Vert\hat{R}\Vert _{\infty}$ and approximate the denominator by $\Vert \hat{R}^{-1} \Vert^{-1}_{\infty}$ in order to compute SERRBD; STRCON returns an approximation RCOND to $1/ (\Vert \hat{R}^{-1} \Vert _{\infty} \Vert \hat{R} \Vert _{\infty})$.

We assume that the rank r of G equals n, because otherwise the $\alpha_i$s and $\beta_i$s are not well determined. For example, if

$$A = \left( \begin{array}{cc} 10^{-16} & 0 \\ 0 & 1 \end{arra... ...( \begin{array}{cc} 10^{-16} & 0 \\ 0 & 1 \end{array} \right)$$

then A and B have $\alpha_{1,2} = 1, .7071$ and $\beta_{1,2} = 0, .7071$, whereas A' and B' have $\alpha'_{1,2} = 0, .7071$ and $\beta'_{1,2} = 1, .7071$, which are completely different, even though $\Vert A - A' \Vert = 10^{-16}$ and $\Vert B - B' \Vert = 10^{-16}$. In this case, $\sigma_{\min} (G) = 10^{-16}$, so G is nearly rank-deficient.

The reason the code fragment assumes $m \geq n$ is that in this case $\hat{R}$ is stored overwritten on A, and can be passed to STRCON in order to compute RCOND. If $m \leq n$, then the first m rows of $\hat{R}$ are stored in A, and the last n-m rows of $\hat{R}$ are stored in B. This complicates the computation of RCOND: either $\hat{R}$ must be copied to a single array before calling STRCON, or else the lower level subroutine SLACON must be used with code capable of solving linear equations with $\hat{R}$ and $\hat{R}^T$ as coefficient matrices.