Reorthogonalization for the Golub–Kahan–Lanczos bidiagonal reduction

Abstract

The Golub–Kahan–Lanczos (GKL) bidiagonal reduction generates, by recurrence, the matrix factorization of $X \in \mathbb{R }^{m \times n}, m \ge n$, given by

$$\begin{aligned} X = UBV^T \end{aligned}$$

where $U \in \mathbb{R }^{m \times n}$ is left orthogonal, $V \in \mathbb{R }^{n \times n}$ is orthogonal, and $B \in \mathbb{R }^{n \times n}$ is bidiagonal. When the GKL recurrence is implemented in finite precision arithmetic, the columns of $U$ and $V$ tend to lose orthogonality, making a reorthogonalization strategy necessary to preserve convergence of the singular values. The use of an approach started by Simon and Zha (SIAM J Sci Stat Comput, 21:2257–2274, 2000) that reorthogonalizes only one of the two left orthogonal matrices $U$ and $V$ is shown to be very effective by the results presented here. Supposing that $V$ is the matrix reorthogonalized, the reorthogonalized GKL algorithm proposed here is modeled as the Householder Q–R factorization of $\left( \begin{array}{c} 0_{n \times k} \\ X V_k \end{array}\right)$ where $V_k = V(:,1:k)$. That model is used to show that if $\varepsilon _M$ is the machine unit and

$$\begin{aligned} \bar{\eta }= \Vert \mathbf{tril }(I-V^T\!~V)\Vert _F, \end{aligned}$$

where $\mathbf{tril }(\cdot )$ is the strictly lower triangular part of the contents, then: (1) the GKL recurrence produces Krylov spaces generated by a nearby matrix $X + \delta X$, $\Vert \delta X\Vert _F = \mathcal O (\varepsilon _M + \bar{\eta }) \Vert X\Vert _F$; (2) singular values converge in the Lanczos process at the rate expected from the GKL algorithm in exact arithmetic on a nearby matrix; (3) a new proposed algorithm for recovering leading left singular vectors produces better bounds on loss of orthogonality and residual errors.

Appendix: Proof of bound on \(\Vert \mathbf{r }_k\Vert _2\)

Proof of Lemma 3.1

We have that, in floating point arithmetic,

$$\begin{aligned} \mathbf{r }_k&= \mathbf{fl }(X^T \mathbf{u }_k -\gamma _k \mathbf{v }_k )\\&= X^T \mathbf{u }_k -\gamma _k \mathbf{v }_k + \delta \mathbf{r }_k \end{aligned}$$

where

$$\begin{aligned} \Vert \delta \mathbf{r }_k\Vert _2&\le \varepsilon _M ((m+1) \Vert X\Vert _F\Vert \mathbf{u }_k\Vert _2 + |\gamma _k|\Vert \mathbf{v }_k\Vert _2 ) {+ \mathcal O (\varepsilon _M ^2 )}\\&= \varepsilon _M ((m+1) \Vert X\Vert _F + |\gamma _k|) {+ \mathcal O (\varepsilon _M ^2 )}. \end{aligned}$$

We note that from (4.39) and (4.40), we have

$$\begin{aligned} \overline{W}_k \left( \begin{array}{c} B_k \\ 0 \end{array}\right) \mathbf{e }_k = \left( \begin{array}{c} 0 \\ X\mathbf{v }_k \end{array}\right) + \delta C_k \mathbf{e }_k \end{aligned}$$

and thus (using the convention that $\phi _1 =0$),

$$\begin{aligned} |\gamma _k |&\le \left\| \left( \begin{array}{c} \phi _k \\ \gamma _k \end{array}\right) \right\| _{2}\\&\le \Vert X\mathbf{v }_k\Vert _2 +\Vert \delta \overline{C}_k \mathbf{e }_k\Vert _2\\&\le \Vert X \mathbf{v }_k\Vert _2 + \Vert X\Vert _F [ f_3(m,n,k)\varepsilon _M + f_2(k) \bar{\eta }_k] + \mathcal O (\varepsilon _M ^2 +\bar{\eta }_k^2) \\&\le \Vert X \Vert _2 + \Vert X\Vert _F [ f_3(m,n,k)\varepsilon _M + f_2(k) \bar{\eta }_k] + \mathcal O (\varepsilon _M ^2 +\bar{\eta }_k^2). \end{aligned}$$

Therefore, since $\varepsilon _M \bar{\eta }_k = \mathcal O (\varepsilon _M ^2 + \bar{\eta }_k^2)$

$$\begin{aligned} \Vert \delta \mathbf{r }_k\Vert _2 \le \varepsilon _M (m+2) \Vert X\Vert _F +\mathcal O (\varepsilon _M ^2+ \bar{\eta }_k^2). \end{aligned}$$

Thus

$$\begin{aligned} \Vert \mathbf{r }_k \Vert _2&\le \Vert X^T \mathbf{u }_k\Vert _2 + |\gamma _k| + \Vert \delta \mathbf{r }_k\Vert _2 \\&\le 2\Vert X\Vert _2+(\varepsilon _M (m+2)+ f_3(m,n,k)\varepsilon _M + f_2(k) \bar{\eta }_k )\Vert X\Vert _F+ \mathcal O (\varepsilon _M ^2+ \bar{\eta }_k^2)) \\&= 2\Vert X\Vert _2 + \mathcal O (\varepsilon _M +\bar{\eta }_k).\\ \end{aligned}$$

$\square$