Exact-arithmetic LSQR with Golub–Kahan bidiagonalization and explicit statuses

Let $A\in {\mathbb{R}}^{m\times n}$, $b\in {\mathbb{R}}^m$, $x_0\in {\mathbb{R}}^n$, $\delta \ge 0$, $\tau >1$, and a budget $N\in \mathbb{N}$ be given. All operations are exact. In the projected matrices and vectors, use local one-based labels: entry $(i,k)$ means the underlying zero-based entry $(i-1,k-1)$, and vector entry $i$ means entry $i-1$. Put $r_0=b-A{x}_0$. If $\mathrm{\parallel }{r}_0{\mathrm{\parallel }}_2=0$, return $(x_0,0,\text{zero-residual})$. If $\mathrm{\parallel }{r}_0{\mathrm{\parallel }}_2\le \tau \delta$, return $(x_0,0,\text{discrepancy})$. If $N=0$, return $(x_0,0,\text{budget-exhausted})$. Otherwise initialize ${\beta }_1=\mathrm{\parallel }{r}_0{\mathrm{\parallel }}_2$, $u_1=r_0\mathrm{/}{\beta }_1$, and ${\alpha }_1=\mathrm{\parallel }{A}^T{u}_1{\mathrm{\parallel }}_2$. If ${\alpha }_1=0$, return $(x_0,0,\text{stationary-nonzero-residual})$. Otherwise set $v_1=A^T{u}_1\mathrm{/}{\alpha }_1$.

Initialize the current best iterate and index to $(x_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}},j_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}})=(x_0,0)$. If a Givens rotation is undefined, return $(x_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}},j_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}},\text{small-QR-breakdown})$.

For $j=1,\dots ,N$, perform these steps.

1. Compute $p_j=A{v}_j-{\alpha }_j{u}_j$ and ${\beta }_{j+1}=\mathrm{\parallel }{p}_j{\mathrm{\parallel }}_2$. If ${\beta }_{j+1}>0$, set $u_{j+1}=p_j\mathrm{/}{\beta }_{j+1}$.
2. If ${\beta }_{j+1}>0$, compute $q_j=A^T{u}_{j+1}-{\beta }_{j+1}{v}_j$ and ${\alpha }_{j+1}=\mathrm{\parallel }{q}_j{\mathrm{\parallel }}_2$. If ${\alpha }_{j+1}>0$, set $v_{j+1}=q_j\mathrm{/}{\alpha }_{j+1}$. A zero ${\beta }_{j+1}$ or ${\alpha }_{j+1}$ marks bidiagonalization breakdown.
3. Put $V_j=[v_1,\dots ,v_j]$. Form the $(j+1)\times j$ lower-bidiagonal matrix $B_j$ with diagonal ${\alpha }_1,\dots ,{\alpha }_j$ and subdiagonal ${\beta }_2,\dots ,{\beta }_{j+1}$, using ${\beta }_{j+1}=0$ if step 1 broke down. Let $e_1=(1,0,\dots ,0{)}^T\in {\mathbb{R}}^{j+1}$ and $c={\beta }_1{e}_1$. For $i=1,\dots ,j$, let $a,b$ be the current entries at positions $(i,i),(i+1,i)$ of the rotated matrix and put $\rho =(a^2+b^2{)}^{1\mathrm{/}2}$. If $\rho =0$, return the current best output with status small-QR-breakdown. Otherwise put $c_i=a\mathrm{/}\rho$, $s_i=b\mathrm{/}\rho$ and apply $\left(\begin{array}{cc}{c}_i& s_i\\ -s_i& c_i\end{array}\right)$ to rows $i,i+1$ of the full matrix and right-hand side. After these $j$ rotations, let $R_j$ be the first $j$ rows of the rotated matrix and $d_j$ the first $j$ entries of the rotated right-hand side. Solve $R_j{z}_j=d_j$ by back substitution. Update $(x_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}},j_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}})=(x_0+V_j{z}_j,j)$.
4. Put $x_j=x_0+V_j{z}_j$. Compute the physical residual $\mathrm{\parallel }A{x}_j-b{\mathrm{\parallel }}_2$. If it is at most $\tau \delta$, return $(x_j,j,\text{discrepancy})$. If bidiagonalization broke down, return $(x_j,j,\text{bidiagonalization-breakdown})$. If $j=N$, return $(x_j,j,\text{budget-exhausted})$; otherwise continue.

The procedure makes no claim that a finite-budget status is a successful regularized approximation. Its exact-arithmetic residual-minimization property is a separate theorem.