# Estimating the error in CG-like algorithms
for least-squares and least-norm problems

Jan Papež  Institute of Mathematics of the Czech Academy of Sciences, Prague,
Czech Republic  Petr Tichý  Faculty of Mathematics and Physics, Charles
University, Prague, Czech Republic

(version of )

###### Abstract

In [Meurant, Papež, Tichý; Numerical Algorithms 88, 2021], we presented an
adaptive estimate for the energy norm of the error in the conjugate gradient
(CG) method. In this paper, we extend the estimate to algorithms for solving
linear approximation problems with a general, possibly rectangular matrix that
are based on applying CG to a system with a positive (semi-)definite matrix
build from the original matrix. We show that the resulting estimate preserves
its key properties: it can be very cheaply evaluated, and it is numerically
reliable in finite-precision arithmetic under some mild assumptions. We
discuss algorithms based on Hestenes–Stiefel-like implementation (often called
CGLS and CGNE in the literature) as well as on bidiagonalization (LSQR and
CRAIG), and both unpreconditioned and preconditioned variants. The numerical
experiments confirm the robustness and very satisfactory behaviour of the
estimate.

## Introduction

Solving linear approximation problems (with a general, possibly rectangular
matrix) is a common task in scientific computing. In this paper, we consider a
least-squares problem

$\min_{z\in\mathbb{R}^{n}}\|b-Az\|,\qquad A\in\mathbb{R}^{m\times n},\
b\in\mathbb{R}^{m},\ m\geq n,$

and a least-norm problem

$\min_{z\in\mathbb{R}^{n}}\|z\|\quad\mbox{subject to}\quad Az=b,\qquad
A\in\mathbb{R}^{m\times n},\ b\in\mathbb{R}^{m},\ b\in\mathcal{R}(A),$

where $\|\cdot\|$ denotes the Euclidean norm and $\mathcal{R}(A)$ is the
_range_ of $A$, i.e., the space spanned by columns of $A$.

When $A$ is large and sparse, an iterative method is typically used to get an
approximation to a (possibly non-unique) solution of the problem. For using
any iterative method, one needs an appropriate stopping criterion which
typically requires an estimate of some quantity that measures the quality of
the computed approximation.

Krylov subspace methods are among the most widely used and studied methods for
iterative solution of least-squares and least-norm problems. Their
mathematical derivation is based on a minimization of a certain error or
residual norm over a Krylov subspace (of a successively increasing dimension);
see, e.g, [21]. Typically, there are several algorithms that are
mathematically equivalent (they are based on the same minimization problem
over the same subspaces) but they differ in the implementation, potentially
leading to a different behavior in finite-precision arithmetic.

The purpose of this paper is to extend a (heuristic-based) adaptive error
estimate derived in [24] for the conjugate gradient (CG) method to a class of
algorithms for solving the above problems. We consider various algorithmic
variants that are mathematically equivalent to applying CG to a system with a
symmetric positive (semi-)definite matrix $A^{T}A$ or $AA^{T}$ build from the
original matrix $A$. In particular, the algorithms discussed in this paper are
CGLS, LSQR, CGNE, and CRAIG (Craig’s method based on bidiagonalization). For
each of these algorithms we present the derivation of the estimate and show
that the two key properties of the estimate known from CG persist. First, the
estimate can be evaluated very cheaply. Second, under the assumption that the
local orthogonality is preserved during finite-precision computations, and
that the maximal attainable accuracy has not yet been reached, the estimate is
numerically reliable even if the convergence of the algorithms is strongly
influenced by rounding errors. The convergence may be significantly delayed
but the delay is reflected in the estimate; it estimates the actual error of
the computed approximation. The rounding error analysis concerning local
orthogonality for the CG algorithm can be found in [32]; see also [23]. Such
an analysis would be doable also for the algorithms discussed in this paper,
however, it would require many additional technical details. The maximal
attainable accuracy was analyzed, e.g., in [16] and [31] in a more general
context. Note that some of the estimates discussed in this paper can be found
in the literature, see, e.g., [19, 3, 2, 11], but without using our adaptive
technique and also without any discussion on using the estimates in finite-
precision arithmetic.

The paper is organized as follows. We first briefly recall the adaptive error
estimate of [24]. Then, in Section 2 we present the least-squares problem and
its solution by CGLS and LSQR, proposed in [18] and [27], respectively. In
Section 3, we discuss the least-norm problem and the error estimate for
Craig’s method (proposed originally in [8, 9], implemented using CGNE ([14, p.
504]), and the algorithm based on Golub–Kahan bidiagonalization [27] that we
call CRAIG; see, e.g., [30]. Then, in Section 4, we discuss the estimate for
preconditioned variants of the algorithms. The results of numerical
experiments are given in Section 5 and the paper ends with concluding remarks.
The MATLAB codes of the algorithms with error estimates are available from the
GitHub repository [28].

##  1 Adaptive error estimate in CG

In this section we briefly recall the adaptive error estimate derived in [24]
for the conjugate gradient method.

First recall the idea of the CG method for solving a system $Ax=b$ with a
symmetric and positive definite matrix $A$. Starting with an initial guess
$x_{0}$ and the associated residual $r_{0}=b-A{}x_{0}$, CG generates the
approximations

$x_{k}\in
x_{0}+\mathcal{K}_{k}(A,r_{0}),\qquad\mathcal{K}_{k}(A,r_{0})\equiv\mbox{span}\\{r_{0},Ar_{0},\ldots,A^{k-1}r_{0}\\},$

such that

$x_{k}\,=\operatorname*{arg\,min}_{y\in
x_{0}+\mathcal{K}_{k}(A,r_{0})}\|x-y\|_{A},$ (1)

where $\|\cdot\|_{A}$ denotes the $A$-norm, $\|v\|_{A}^{2}=v^{T}Av$. The
standard Hestenes and Stiefel implementation of the CG method is given in
Algorithm 1. Therein, we put on line 6 a call of the function
$\mathtt{adaptive}(\cdot)$ that performs the computation of the adaptive error
estimate and is explained below.

Note that CG can be applied also to a system $Ax=b$ with a matrix $A$ that is
symmetric positive semi-definite only; see, e.g., [20]. Then $\|\cdot\|_{A}$
is a norm on $\mathcal{R}(A)$ only, and we need an additional assumption so
that (1) makes sense. In particular, we must assume that the system is
consistent, i.e., that $b\in\mathcal{R}(A)$, and choose the initial guess
$x_{0}$ from $\mathcal{R}(A)$. This assures that
$\mathcal{K}_{k}(A,r_{0})\subset\mathcal{R}(A)$ and $x_{k}\in\mathcal{R}(A)$,
for $k=0,1,\ldots$, so that the minimization (1) can be considered over
$\mathcal{R}(A)$. Note that in finite-precision arithmetic, the singularity of
the system matrix $A$ can cause some numerical instabilities, in particular,
the error norm can start to increase after reaching the level of maximal
attainable accuracy; see [20, Section 4].

Algorithm 1 Conjugate Gradients for $Ax=b$ (CG)

1: input $A$, $b$, $x_{0}$

2: $r_{0}=b-Ax_{0}$, $p_{0}=r_{0}$

3: $\ell=0$, $\tau=0.25$prescribe the tolerance for error estimation

4: for $k=0,1,\dots$ until convergence do

5:
$\gamma_{k}={\|r_{k}\|^{2}}/{\|p_{k}\|_{A}^{2}}={r_{k}^{T}r_{k}}\,/\,{p_{k}^{T}{A}p_{k}}$

6: ${\tt adaptive}\/(k,\ell,\Delta_{k}=\gamma_{k}\|r_{k}\|^{2},\tau)$

7: $x_{k+1}=x_{k}+\gamma_{k}p_{k}$

8: $r_{k+1}=r_{k}-\gamma_{k}Ap_{k}$

9: $\delta_{k+1}=\|r_{k+1}\|^{2}/\|r_{k}\|^{2}$

10: $p_{k+1}=r_{k+1}+\delta_{k+1}p_{k}$

11: end for

The adaptive error estimate of [24] is based on the formula (which we call the
Hestenes–Stiefel formula),

$\|x-x_{\ell}\|^{2}_{A}=\sum^{k}_{j=\ell}\Delta_{j}+\|x-x_{k+1}\|^{2}_{A},\qquad\Delta_{j}\equiv\gamma_{j}\|r_{j}\|^{2},$
(2)

for $\ell\leq k$, and on a procedure that adaptively finds $\ell=\ell(k)$, for
a given $k$, such that

$\frac{\|x-x_{\ell}\|^{2}_{A}-\sum^{k}_{j=\ell}\Delta_{j}}{\|x-x_{\ell}\|^{2}_{A}}\leq\tau$

for some prescribed tolerance $\tau\in(0,1)$ (typically set as $0.25$). Then
the error $\|x-x_{k}\|^{2}_{A}$ is estimated by

$\Delta_{\ell:k}\equiv\sum^{k}_{j=\ell}\Delta_{j}\approx\|x-x_{\ell}\|^{2}_{A}.$

From (2), it is clear that $\Delta_{\ell:k}$ represents a lower bound on
$\|x-x_{\ell}\|^{2}_{A}$. The bound $\Delta_{\ell:k}$ is very cheap to
evaluate, it is numerically stable [32, 33], and if $\ell$ is set properly, it
provides a sufficiently accurate estimate.

Algorithm 2 ${\mathtt{adaptive}()}$ — adaptive choice of the delay $k-\ell$

1: input $k$, $\ell$, $\Delta_{k}$, $\tau$

2: $\mathrm{TOL}=10^{-4}$set a parameter for error estimation

3: output $\ell$

4: set $m$ as the largest index $j$, $0\leq j<k$, such that
$\frac{\Delta_{\ell:k}}{\Delta_{j:k}}\leq\mbox{TOL}$ if such an $m$ does not
exist, set $m=0$

5: determine $S$ as  $S=\max_{m\leq j<k}\Delta_{j:k}/\Delta_{j}$

6: while $S\Delta_{k}/\Delta_{\ell:k-1}\leq\tau$ and $\ell<k$ do

7: accept $\Delta_{\ell:k}$ as an estimate of the error at $\ell$th iteration

8: $\ell=\ell+1$

9: end while

The function $\mathtt{adaptive}(\cdot)$ for setting the (adaptive) delay
$k-\ell$ and computing the error estimate is presented in Algorithm 2. As one
can observe, in the $k$th CG iteration when the approximation $x_{k+1}$ is
computed, the error estimate is evaluated for some previous approximation
$x_{\ell}$. For a detailed description, derivation, and reasoning, see [24].
Note that in the original paper [24], $\Delta_{\ell:k-1}$ was used on line 7
of Algorithm 2 for estimating the error. Since $\Delta_{\ell:k}$ is available
and $\Delta_{\ell:k-1}<\Delta_{\ell:k}\leq\|x-x_{\ell}\|^{2}_{A}$, we suggest
using a tighter estimate. Obviously, this function can be easily added to the
existing CG codes.

As described in [24, Section 4], the adaptive error estimate can also be
extended to preconditioned CG while preserving all key properties. For the
sake of simplicity, the resulting algorithm is not presented here and we refer
to [24, Section 4].

##  2 Estimating the error in least-squares problems

Let us consider a least-squares problem

$x=\operatorname*{arg\,min}_{z\in\mathbb{R}^{n}}\|b-Az\|$ (3)

where $A\in\mathbb{R}^{m\times n}$, $m\geq n$, is a given matrix and
$b\in\mathbb{R}^{m}$ is a right-hand side. It is well known that $x$ is a
least-squares solution if and only if $b-Ax\perp\mathcal{R}(A)$, or
equivalently, if and only if $x$ is the solution of the system of _normal
equations_

$A^{T}Ax=A^{T}b.$ (4)

If $\mathrm{rank}(A)=n$, then the matrix $A^{T}A$ is nonsingular and the
solution of (3) is unique. Moreover, $A^{T}A$ is symmetric and positive
definite and, therefore, the CG method can be applied to (4).

When CG is applied to (4) starting with an initial guess $x_{0}$ and
$r_{0}=b-Ax_{0}$, the approximation $x_{k}$ computed at the $k$th iteration
satisfies

$x_{k}\in x_{0}+\mathcal{K}_{k}(A^{T}A,A^{T}r_{0}),$

and

$x_{k}=\operatorname*{arg\,min}_{y\in
x_{0}+\mathcal{K}_{k}(A^{T}A,A^{T}r_{0})}\|x-y\|_{A^{T}A}.$ (5)

There are several ways how CG for (4) can be implemented. For an overview we
refer to [5] or to a series of papers [11, 12, 13].

Algorithm 3 CGLS

1: input $A$, $b$, $x_{0}$

2: $r_{0}=b-Ax_{0}$

3: $s_{0}=p_{0}=A^{T}r_{0}$

4: for $k=0,1,\dots$ do

5: $q_{k}=Ap_{k}$

6: $\gamma_{k}=\|s_{k}\|^{2}/\|q_{k}\|^{2}$

7: $x_{k+1}=x_{k}+\gamma_{k}p_{k}$

8: $r_{k+1}=r_{k}-\gamma_{k}q_{k}$

9: $s_{k+1}=A^{T}r_{k+1}$

10: $\delta_{k+1}=\|s_{k+1}\|^{2}/\|s_{k}\|^{2}$

11: $p_{k+1}=s_{k+1}+\delta_{k+1}p_{k}$

12: end for

The CGLS algorithm (Algorithm 3) has been described already in the CG seminal
paper [18, p. 424]. It corresponds to applying the CG algorithm to the system
of normal equations with a slight algebraic rearrangement to avoid the vectors
of the form $A^{T}Ap$. Due to this rearrangement, CGLS has better numerical
properties than CG (Algorithm 1) naively applied to (4). The CGLS algorithm
can be found, e.g., in [27, Section 7.1] and [5, Algorithm 3.1].

Algorithm 4 LSQR

1: input $A$, $b$

2: $\beta_{1}u_{1}=b$$\beta_{i}$ are the normalization coefficients to have
$\|u_{i}\|=1$, $i=1,2,\ldots$

3: $\alpha_{1}v_{1}=A^{T}u_{1}$ $\alpha_{i}$ are the normalization
coefficients to have $\|v_{i}\|=1$, $i=1,2,\ldots$

4: $w_{1}=v_{1}$

5: $x_{0}$ = 0

6: $\bar{\phi}_{1}=\beta_{1}$

7: $\bar{\rho}_{1}=\alpha_{1}$

8: for $k=1,2,\dots$ do

9: $\beta_{k+1}u_{k+1}=Av_{k}-\alpha_{k}u_{k}$

10: $\alpha_{k+1}v_{k+1}=A^{T}u_{k+1}-\beta_{k+1}v_{k}$

11: $\rho_{k}=(\bar{\rho}_{k}^{2}+\beta^{2}_{k+1})^{1/2}$

12: $c_{k}=\bar{\rho}_{k}/\rho_{k}$

13: $s_{k}=\beta_{k+1}/\rho_{k}$

14: $\theta_{k+1}=s_{k}\alpha_{k+1}$

15: $\bar{\rho}_{k+1}=-c_{k}\alpha_{k+1}$

16: $\phi_{k}=c_{k}\bar{\phi}_{k}$

17: $\bar{\phi}_{k+1}=s_{k}\bar{\phi}_{k}$

18: $x_{k}=x_{k-1}+(\phi_{k}/\rho_{k})w_{k}$

19: $w_{k+1}=v_{k+1}+(\theta_{k+1}/\rho_{k})w_{k}$

20: end for

A mathematically equivalent algorithm to CGLS based on Golub–Kahan
bidiagonalization has been derived in [27]. The algorithm is called LSQR and
it is given in Algorithm 4.

Since CG can also be applied to a singular system (see the corresponding
discussion in Section 1), the CGLS and LSQR algorithms can also be applied if
$\mathrm{rank}(A)<n$. In such case, however, the least-squares solution of (3)
is not unique.

###  2.1 Error estimation in CGLS and LSQR

Let $b_{|\mathcal{R}(A)}$ be the orthogonal projection of $b$ onto the range
of $A$, and denote

$r=b-b_{|\mathcal{R}(A)}.$

In this section we focus on the estimation of the error

$\|x-x_{k}\|^{2}_{A^{T}A}=\|A(x-x_{k})\|^{2}=\|b-r-
Ax_{k}\|^{2}=\|r_{k}-r\|^{2}=\|r_{k}\|^{2}-\|r\|^{2},$

where we set $r_{k}=b-Ax_{k}$ and used the fact that $r\perp r-r_{k}$.
Therefore, if $b\in\mathcal{R}(A)$, which means that $r=0$, the error
$\|x-x_{k}\|^{2}_{A^{T}A}$ is equal to the fully computable quantity
$\|r_{k}\|^{2}=\|b-Ax_{k}\|^{2}$ and there is no need for an estimator.
Hereafter, we therefore assume that $b\notin\mathcal{R}(A)$. Then the norm
$\|r\|$ is non-zero and unknown, and the residual norm $\|r_{k}\|$ alone may
not provide a sufficient information to set a proper stopping criterion. This
has been discussed, e.g., in [19]. It seems natural to set a stopping
criterion based on comparison of (the norms of) the residual $r$ and the
iterative residual $r_{k}$. Such discussion is, however, beyond the scope of
this paper.

The error norm $\|x-x_{k}\|_{A^{T}A}$ is the relevant quantity for the CGLS
and LSQR algorithms as it is minimized within the iterations (it is therefore
monotone) and, as we will see in Section 4, this norm is also minimized if a
preconditioning is considered. In many practical situations, the users are
also interested in estimating the Euclidean norm of the error $\|x-x_{k}\|$.
This error norm can be efficiently estimated in CG without preconditioning,
based on a direct relation between the Euclidean and the energy norm of the
error, see, e.g. [18] or [32], or, based on the relation of CG with Gauss
quadrature [15, 22]. Then similar techniques for improving the accuracy of the
estimate, that use an analog of (2), can be developed. All these results are
transferable from CG to CGLS and LSQR, but, again, without preconditioning.
When using preconditioning, there is no direct relation between the Euclidean
and the energy norm of the error in general. However, one can bound the
Euclidean norm of the error using

$\|x-x_{k}\|\leq\frac{1}{\sigma}\,\|x-x_{k}\|_{A^{T}A}$

that requires an a priori knowledge of a lower bound $\sigma$ on the smallest
singular value of $A$; see, e.g., [2, 13, 17]. Clearly, the above upper bound
need not represent a good estimate of $\|x-x_{k}\|$, but, in general, we do
not have anything better in hands. Note that to use the above upper bound, we
still need to have an estimate of the quantity $\|x-x_{k}\|_{A^{T}A}$, which
is our aim in this paper.

###  2.2 Error estimate in CGLS

The first step to derive an estimate as in [24] is to find an expression
analogous to (2). By a simple algebraic manipulation, see [33, p. 798], we
obtain

$\left\|x-x_{k}\right\|_{A^{T}A}^{2}-\left\|x-x_{k+1}\right\|_{A^{T}A}^{2}=\|x_{k+1}-x_{k}\|_{A^{T}A}^{2}+2\left(x_{k+1}-x_{k}\right)^{T}A^{T}A\left(x-x_{k+1}\right).$
(6)

In the following we use only relations that do hold (up to a small inaccuracy)
also during computations in finite-precision arithmetic. To shorten the
terminology, we say that these identities hold numerically. In particular, we
avoid using global orthogonality of vectors that is usually lost quickly.
Since in CGLS

$x_{k+1}-x_{k}=\gamma_{k}p_{k},\qquad s_{k+1}=A^{T}A(x-x_{k+1})$ (7)

we obtain

$\left\|x-x_{k}\right\|_{A^{T}A}^{2}-\left\|x-x_{k+1}\right\|_{A^{T}A}^{2}=\gamma_{k}^{2}\|p_{k}\|_{A^{T}A}^{2}+2\gamma_{k}p_{k}^{T}s_{k+1}.$

Using the formula for computing $\gamma_{k}$ (line 6) and the (local)
orthogonality between $s_{k+1}$ and $p_{k}$,

$\gamma_{k}=\frac{\|s_{k}\|^{2}}{\|p_{k}\|_{A^{T}A}^{2}},\quad
p_{k}^{T}s_{k+1}=0,$ (8)

we get an analog to (2),

$\|x-x_{\ell}\|^{2}_{A^{T}A}=\sum_{j=\ell}^{k}\gamma_{j}\|s_{j}\|^{2}+\|x-x_{k+1}\|^{2}_{A^{T}A}$

and the error estimator

$\Delta^{\mathrm{CGLS}}_{\ell:k}\equiv\sum_{j=\ell}^{k}\gamma_{j}\|s_{j}\|^{2}\approx\|x-x_{\ell}\|^{2}_{A^{T}A}=\|r_{\ell}\|^{2}-\|r\|^{2}.$

To derive the analog of (2) we have used relations (7) and (8), which do hold
numerically until the level of maximal attainable accuracy is reached; see,
e.g., [32].

###  2.3 Error estimate in LSQR

For LSQR, the situation is more delicate. Here the explanation needs more
space and goes hand in hand with a derivation, inspired by ideas of [32].

The vectors $u_{i}$ and $v_{i}$ generated in LSQR (Algorithm 4) by Golub–Kahan
bidiagonalization satisfy

$\displaystyle AV_{k}$ $\displaystyle=$ $\displaystyle U_{k+1}B_{k}$
$\displaystyle A^{T}U_{k+1}$ $\displaystyle=$ $\displaystyle
V_{k}B_{k}^{T}+\alpha_{k+1}v_{k+1}e_{k+1}^{T}$

where $V_{k}=[v_{1},v_{2},\ldots,v_{k}]$,
$U_{k+1}=[u_{1},u_{2},\ldots,u_{k+1}]$, and

$B_{k}=\left[\begin{array}[]{ccccc}\alpha_{1}\\\ \beta_{2}&\ddots\\\
&\ddots&\ddots\\\ &&\beta_{k-1}&\alpha_{k-1}\\\ &&&\beta_{k}&\alpha_{k}\\\
&&&&\beta_{k+1}\end{array}\right].$

By a simple algebraic manipulation we get

$A^{T}AV_{k}=V_{k}T_{k}+\alpha_{k+1}\beta_{k+1}v_{k+1}e_{k+1}^{T},\qquad
T_{k}=B_{k}^{T}B_{k}.$

so that $v_{i}$ can be seen as Lanczos vectors generated by the corresponding
three-term recurrence for $A^{T}A$.

Assuming $x_{0}=0$, the LSQR approximation is given by

$x_{k}=V_{k}y_{k},\quad y_{k}=\arg\min\|\beta_{1}e_{1}-B_{k}y\|$ (9)

and the corresponding least-squares problem is solved using the QR
factorization of $B_{k}$. In particular,

$Q_{k}\begin{bmatrix}B_{k}&\beta_{1}e_{1}\end{bmatrix}=\begin{bmatrix}R_{k}&f_{k}\\\
&\tilde{\phi}_{k+1}\end{bmatrix}=\left[\begin{array}[]{cccccc}\rho_{1}&\theta_{2}&&&&\phi_{1}\\\
&\rho_{2}&\theta_{3}&&&\phi_{2}\\\ &&\ddots&\ddots&&\vdots\\\
&&&\rho_{k-1}&\theta_{k}&\phi_{k-1}\\\ &&&&\rho_{k}&\phi_{k}\\\
&&&&&\tilde{\phi}_{k+1}\end{array}\right],$

and $y_{k}$ is the solution of $R_{k}y=f_{k}$. Therefore,
$T_{k}=B_{k}^{T}B_{k}=R_{k}^{T}R_{k}$, so that $R_{k}^{T}R_{k}$ represents
Cholesky factorization of $T_{k}$. Since CGLS (CG applied to $A^{T}Ax=A^{T}b$)
computes (implicitly) also the Cholesky factorization of $T_{k}$, see, e.g.,
[25], we get the relation among coefficients that appear in both algorithms,

$\left[\begin{array}[]{cccc}\frac{1}{\sqrt{\gamma_{0}}}\\\
\sqrt{\frac{\delta_{1}}{\gamma_{0}}}&\ddots\\\ &\ddots&\ddots\\\
&&\sqrt{\frac{\delta_{k-1}}{\gamma_{k-2}}}&\frac{1}{\sqrt{\gamma_{k-1}}}\end{array}\right]=\left[\begin{array}[]{ccccc}\rho_{1}\\\
\theta_{2}&\ddots\\\ &&\ddots\\\ &&\theta_{k-1}&\rho_{k-1}\\\
&&&\theta_{k}&\rho_{k}\end{array}\right]$ (10)

where $\gamma_{i}$ and $\delta_{i}$ are the CGLS coefficients; see Algorithm
3.

It is well known that during finite-precision computations, the global
orthogonality among the vectors in columns of $V_{k}$ and $U_{k}$ is usually
lost very quickly. In other words, one cannot expect that
$V_{k}^{T}V_{k}=I=U_{k}^{T}U_{k}$. Our aim is to explain that despite the loss
of global orthogonality and under some natural assumptions which will be
specified later, the Hestenes–Stiefel formula for LSQR, which has the form

$\left\|x-x_{k-1}\right\|_{A^{T}A}^{2}-\left\|x-x_{k}\right\|_{A^{T}A}^{2}=\phi_{k}^{2},$
(11)

still holds numerically. Note that we cannot directly apply results of [32,
33] since LSQR (Algorithm 4) uses different recurrences than CG (Algorithm 1).
We do not present here a detailed rounding error analysis for the formula (11)
like in [32, 33]. Instead, we focus on a justification based on a convenient
derivation of the formula (11) and the discussion about preserving local
orthogonality. We start with the following lemma; see also [27, p. 52].

######  Lemma 1.

For the quantities generated by LSQR (Algorithm 4) it holds that

$\displaystyle A^{T}(b-Ax_{k})$ $\displaystyle=$
$\displaystyle-\theta_{k+1}\phi_{k}v_{k+1}.$ (12)

###### Proof.

Note that $y_{k}$ is the solution of the least-squares problem (9) so that
$B_{k}^{T}B_{k}y_{k}=B_{k}^{T}\beta_{1}e_{1}$. Since $x_{k}=V_{k}y_{k}$ and
$AV_{k}=U_{k+1}B_{k}$, we get

$\displaystyle A^{T}(b-Ax_{k})$ $\displaystyle=$ $\displaystyle
A^{T}U_{k+1}\left(\beta_{1}e_{1}-B_{k}y_{k}\right)$ $\displaystyle=$
$\displaystyle\left(V_{k}B_{k}^{T}+\alpha_{k+1}v_{k+1}e_{k+1}^{T}\right)\left(\beta_{1}e_{1}-B_{k}y_{k}\right)$
$\displaystyle=$ $\displaystyle
V_{k}\left(B_{k}^{T}\beta_{1}e_{1}-B_{k}^{T}B_{k}y_{k}\right)-\alpha_{k+1}v_{k+1}e_{k+1}^{T}B_{k}y_{k}$
$\displaystyle=$
$\displaystyle-\alpha_{k+1}v_{k+1}\left[e_{k+1}^{T}B_{k}y_{k}\right]$
$\displaystyle=$
$\displaystyle\alpha_{k+1}v_{k+1}e_{k+1}^{T}\left(\beta_{1}e_{1}-B_{k}y_{k}\right),$

where we have used $e_{k+1}^{T}e_{1}=0$. Using the QR factorization

$\begin{bmatrix}B_{k}&\beta_{1}e_{1}\end{bmatrix}=Q_{k}^{T}\begin{bmatrix}R_{k}&f_{k}\\\
&\bar{\phi}_{k+1}\end{bmatrix}$

and $R_{k}y_{k}=f_{k}$ we obtain

$\beta_{1}e_{1}-B_{k}y_{k}=Q_{k}^{T}\begin{bmatrix}0\\\
\bar{\phi}_{k+1}\end{bmatrix};$

see [27, equation (4.8)]. Hence, the last component of
$\beta_{1}e_{1}-B_{k}y_{k}$ is given by $-c_{k}\bar{\phi}_{k+1}$. Moreover,
using
$\alpha_{k+1}c_{k}\bar{\phi}_{k+1}=\left(\alpha_{k+1}s_{k}\right)\left(c_{k}\bar{\phi}_{k}\right)=\theta_{k+1}\phi_{k}$
we get

$\displaystyle A^{T}(b-Ax_{k})$ $\displaystyle=$
$\displaystyle-\alpha_{k+1}c_{k}\bar{\phi}_{k+1}v_{k+1}=-\theta_{k+1}\phi_{k}v_{k+1}.$

∎

Since the relations $x_{k}=V_{k}y_{k}$ and $AV_{k}=U_{k+1}B_{k}$ hold also for
the computed vectors and coefficients (up to an inaccuracy comparable to
machine precision $\varepsilon$, norms of the vectors, and $\|A\|^{2}$), and
since the computation of $R_{k}$ is backward stable, we can expect that (12)
holds numerically, until the level of maximal attainable accuracy is reached.

Now we derive an analog of (11) without using any orthogonality assumption. We
will see that the resulting formula will contain local orthogonality terms

$\ell_{k}^{(w)}\equiv\frac{\theta_{k+1}}{\rho_{k}}v_{k+1}^{T}w_{k}$ (13)

which are zero in exact arithmetic.

######  Theorem 1.

For the quantities generated by LSQR (Algorithm 4) it holds that

$\left\|x-x_{k-1}\right\|_{A^{T}A}^{2}-\left\|x-x_{k}\right\|_{A^{T}A}^{2}=\phi_{k}^{2}\left({1+\ell_{k-1}^{(w)}}-\ell_{k}^{(w)}\right).$

###### Proof.

By a simple algebraic manipulation we obtain

$\left\|x-x_{k-1}\right\|_{A^{T}A}^{2}-\left\|x-x_{k}\right\|_{A^{T}A}^{2}=\left(x_{k}-x_{k-1}\right)^{T}A^{T}A\left(x-x_{k-1}\right)+\left(x_{k}-x_{k-1}\right)^{T}A^{T}A\left(x-x_{k}\right),$

which is just a minor modification of the formula (6). Using the relations

$A^{T}A(x-x_{k})=-\theta_{k+1}\phi_{k}v_{k+1},\quad
x_{k}-x_{k-1}=\frac{\phi_{k}}{\rho_{k}}w_{k}$

we obtain

$\left\|x-x_{k-1}\right\|_{A^{T}A}^{2}-\left\|x-x_{k}\right\|_{A^{T}A}^{2}=\phi_{k}^{2}\left(\left[-\frac{\theta_{k}\phi_{k-1}}{\rho_{k}\phi_{k}}\right]w_{k}^{T}v_{k}-\frac{\theta_{k+1}}{\rho_{k}}v_{k+1}^{T}w_{k}\right).$

Note that

$-\frac{\theta_{k}\phi_{k-1}}{\rho_{k}\phi_{k}}=-\frac{(s_{k-1}\alpha_{k})(c_{k-1}\bar{\phi}_{k-1})}{\rho_{k}\phi_{k}}=\frac{\bar{\rho}_{k}\bar{\phi}_{k}}{\rho_{k}\phi_{k}}=1$
(14)

so that

$\left\|x-x_{k-1}\right\|_{A^{T}A}^{2}-\left\|x-x_{k}\right\|_{A^{T}A}^{2}=\phi_{k}^{2}\left(w_{k}^{T}v_{k}-\frac{\theta_{k+1}}{\rho_{k}}v_{k+1}^{T}w_{k}\right).$

Finally, from line 19 of Algorithm 4,

$v_{k}^{T}w_{k}=v_{k}^{T}\left(v_{k}+\frac{\theta_{k}}{\rho_{k-1}}w_{k-1}\right)=1+\frac{\theta_{k}}{\rho_{k-1}}v_{k}^{T}w_{k-1}.$

∎

As a consequence, (11) should hold numerically, if the local orthogonality
terms (13) are small (meaning that their magnitude is much less than one) and
if (12) holds numerically (if the level of maximal attainable accuracy have
not been reached yet). Note that (12) has been used in the proof of Theorem 1.
Hence, the problem of justification of the identity (11) in finite-precision
arithmetic is in this way reduced to the problem of bounding local
orthogonality between the computed vectors $v_{k+1}$ and $w_{k}$.

To understand better the terms $\ell_{k}^{(w)}$, let us realize that LSQR
computes in exact arithmetic the same approximations $x_{k}$ as CGLS.
Therefore,

$s_{k}=A^{T}(b-Ax_{k})=-\theta_{k+1}\phi_{k}v_{k+1},$ (15)

where $s_{k}=A^{T}r_{k}$ are computed in Algorithm 3. By comparing
$x_{k}-x_{k-1}$ in both algorithms, we obtain

$\gamma_{k-1}p_{k-1}=\frac{\phi_{k}}{\rho_{k}}w_{k}\,.$

Therefore, CGLS and LSQR compute the same vectors, just scaled differently. In
particular, $w_{k}$ is a multiple of $p_{k-1}$ and $v_{k+1}$ is a multiple of
$s_{k}$. Both algorithms can be seen as a variant of CG applied to the system
of normal equations $A^{T}Ax=A^{T}b$. Finally, realizing that
$\gamma_{k-1}=\rho_{k}^{-2}$, see (10), we obtain

$\ell_{k}^{(w)}=\frac{\theta_{k+1}}{\rho_{k}}v_{k+1}^{T}w_{k}=-\frac{s_{k}^{T}p_{k-1}}{\rho^{2}_{k}\phi^{2}_{k}}=-\frac{s_{k}^{T}p_{k-1}}{\|s_{k-1}\|^{2}},$

where the last equality follows from (14) and (15), which gives
$\|s_{k-1}\|^{2}=\theta^{2}_{k}\phi^{2}_{k-1}$. Therefore, the local
orthogonality term $\ell_{k}^{(w)}$ can be seen as a counterpart of the term
that was analyzed in detail for CG; see [32, Theorem 9.1]. Based on results of
[32] one can expect that the size of $\ell_{k}^{(w)}$ can be bounded by
machine precision $\varepsilon$ multiplied by some factor that can depend on
$\kappa(A^{T}A)$ (or its analog if $A$ is singular), dimension of the problem
and the number of iterations. However, a proper rounding error analysis
leading to the mentioned result is beyond the scope of this paper. Below we
check the validity of (11) numerically, by plotting the relative error

$\frac{\left|\big{(}\left\|x-x_{k-1}\right\|_{A^{T}A}^{2}-\left\|x-x_{k}\right\|_{A^{T}A}^{2}\big{)}-\phi_{k}^{2}\right|}{\left\|x-x_{k-1}\right\|_{A^{T}A}^{2}-\left\|x-x_{k}\right\|_{A^{T}A}^{2}}$
(16)

for the challenging (in terms of numerical stability) test case from [5] with
the conditioning $\kappa(A^{T}A)=10^{12}$; see Figure 1. We observe that the
relative error is (approximately) inversely proportional to the error
$\left\|x-x_{k}\right\|_{A^{T}A}$ and stays significantly below 1 until the
maximal attainable accuracy is reached.

Figure 1: Test case from [5, p. 733]: error norm and the relative error (16).

In summary, (11) results in

$\|x-x_{\ell}\|^{2}_{A^{T}A}=\sum_{j=\ell}^{k}{\phi}_{j+1}^{2}+\|x-x_{k+1}\|^{2}_{A^{T}A}$
(17)

yielding the error estimator

$\Delta^{\mathrm{LSQR}}_{\ell:k}\equiv\sum_{j=\ell}^{k}{\phi}_{j+1}^{2}\approx\|x-x_{\ell}\|^{2}_{A^{T}A}=\|r_{\ell}\|^{2}-\|r\|^{2}.$

The error estimator contains only scalars that are available during the LSQR
iterations and it is very cheap to evaluate. It should be reliable also during
finite-precision computations, if local orthogonality between $v_{k}$ and
$w_{k-1}$ is well preserved and until the level of maximal attainable accuracy
is reached.

##  3 Estimating the error in least-norm problems

Let us now consider a least-norm problem

$x=\operatorname*{arg\,min}_{z\in\mathbb{R}^{n}}\|z\|\quad\mbox{subject
to}\quad Az=b,$ (18)

where $A\in\mathbb{R}^{m\times n}$ and $b\in\mathbb{R}^{m}$,
$b\in\mathcal{R}(A)$. A method of choice for solving this (consistent) problem
might be Craig’s method that is described in the following section.

###  3.1 Craig’s method

The idea of the method is to write the solution $x$ of (18) as $x=A^{T}y$ for
a proper vector $y$ and compute the approximations $x_{k}=A^{T}y_{k}$, where
$y_{k}$ is an approximation to $y$ given by the CG method applied to the
system

$AA^{T}y=b,$ (19)

starting with some $y_{0}$. Denoting $r_{0}=b-AA^{T}y_{0}=b-Ax_{0}$, the
approximations $y_{k}$ satisfy

$y_{k}\in
y_{0}+\mathcal{K}_{k}(AA^{T},r_{0}),\qquad\|y-y_{k}\|_{AA^{T}}=\min_{z\in
y_{0}+\mathcal{K}_{k}(AA^{T},r_{0})}\|y-z\|_{AA^{T}}.$

A simple algebraic manipulation shows that

$\|y-y_{k}\|^{2}_{AA^{T}}=\|x-x_{k}\|^{2},$ (20)

which means that Craig’s method minimizes the Euclidean norm of the error of
$x_{k}$ over the affine space

$A^{T}y_{0}+A^{T}\mathcal{K}_{k}(AA^{T},r_{0})=x_{0}+\mathcal{K}_{k}(A^{T}A,A^{T}r_{0}).$

The above discussion can be applied also if the matrix $AA^{T}$ is singular;
see the discussion in Section 1.

Analogously to solving the normal equations (4), there are several ways how
the CG method for (19) can be applied. A study of numerical stability of
various implementations was presented in [5].

The first algorithm was proposed by Craig in [8, 9], developed originally as
an algorithm for solving unsymmetric systems (see the review in MathSciNet
written by Forsythe). Fadeev and Fadeeva developed almost the same algorithm
[14, p. 504] with the difference that instead of direction vectors, the
algorithm updates the vectors $g_{k}\equiv A^{T}p_{k}$. The resulting
algorithm is called CGNE and is listed in Algorithm 5. The same algorithm, up
to a different notation, is also given in [5, Algorithm 2.1] when setting
$\mu=0$ therein.

Algorithm 5 CGNE

1: input $A$, $b$, $x_{0}$

2: $r_{0}=b-Ax_{0}$

3: $p_{0}=A^{T}r_{0}$

4: for $k=0,1,2,\dots$ do

5: $\gamma{}_{k}={\|r_{k}\|^{2}}/{\|p_{k}\|^{2}}$

6: $x_{k+1}=x_{k}+\gamma{}_{k}p_{k}$

7: $r_{k+1}=r_{k}-\gamma{}_{k}Ap_{k}$

8: $\delta_{k+1}={\|r_{k+1}\|^{2}}/{\|r_{k}\|^{2}}$

9: $p_{k+1}=A^{T}r_{k+1}+\delta_{k+1}p_{k}$

10: end for

Paige [26] and Paige and Saunders [27] developed another mathematically
equivalent version of this algorithm using the Golub–Kahan bidiagonalization.
The algorithm can also compute the approximations $y_{k}$ to $AA^{T}y=b$. We
present it in Algorithm 6 and call it CRAIG, to be consistent with the
notation introduced by Paige and Saunders.

Algorithm 6 CRAIG

1: input $A$, $b$

2: $\zeta_{0}=-1$

3: $x_{0}=y_{0}=d_{0}=0$

4: $\beta_{1}u_{1}=b$ $\beta_{i}$ are the normalization coefficients to have
$\|u_{i}\|=1$, $i=1,2,\ldots$

5: $\alpha_{1}v_{1}=A^{T}u_{1}$ $\alpha_{i}$ are the normalization
coefficients to have $\|v_{i}\|=1$, $i=1,2,\ldots$

6: for $k=1,2,\dots$ do

7: $\zeta_{k}=-\frac{\beta_{k}}{\alpha_{k}}\zeta_{k-1}$

8: $x_{k}=x_{k-1}+\zeta_{k}v_{k}$

9: $d_{k}=\left(u_{k}-\beta_{k}d_{k-1}\right)/\alpha_{k}$

10: $y_{k}=y_{k-1}+\zeta_{k}d_{k}$ lines 9–10 can be removed if the
approximation $y_{k}$ is not needed

11: $\beta_{k+1}u_{k+1}=Av_{k}-\alpha_{k}u_{k}$

12: $\alpha_{k+1}v_{k+1}=A^{T}u_{k+1}-\beta_{k+1}v_{k}$

13: end for

We now present a way to estimate the error $\|x-x_{k}\|^{2}$ in both
implementations of Craig’s method. Similarly to (6), by a simple algebraic
manipulation we obtain

$\left\|x-x_{k}\right\|^{2}-\left\|x-x_{k+1}\right\|^{2}=\|x_{k+1}-x_{k}\|^{2}+2\left(x_{k+1}-x_{k}\right)^{T}\left(x-x_{k+1}\right).$
(21)

###  3.2 Error estimate in CGNE

Let us first realize that, in CGNE,

$\|x_{k+1}-x_{k}\|^{2}=\gamma_{k}^{2}\|p_{k}\|^{2}=\gamma_{k}\|r_{k}\|^{2},$

where we have used $\gamma_{k}\|p_{k}\|^{2}=\|r_{k}\|^{2}$.

Consider formally the recurrence for computing vectors $y_{k}$ such that
$x_{k}=A^{T}y_{k}$, and the recurrence for computing the corresponding
direction vectors $s_{k}$ such that $p_{k}=A^{T}s_{k}$,

$\displaystyle y_{k+1}$ $\displaystyle=$ $\displaystyle
y_{k}+\gamma_{k}s_{k},$ $\displaystyle s_{k+1}$ $\displaystyle=$
$\displaystyle r_{k+1}+\delta_{k+1}s_{k},$

with $y_{0}$ satisfying $x_{0}=A^{T}y_{0}$ and $s_{0}=r_{0}$. Then the last
term on the right hand side of (21) corresponds to

$\left(x_{k+1}-x_{k}\right)^{T}\left(x-x_{k+1}\right)=\left(y_{k+1}-y_{k}\right)^{T}A\left(x-x_{k+1}\right)=\gamma_{k}s_{k}^{T}r_{k+1},$

so that

$\left\|x-x_{k}\right\|^{2}-\left\|x-x_{k+1}\right\|^{2}=\gamma_{k}\|r_{k}\|^{2}+2\gamma_{k}s_{k}^{T}r_{k+1}.$

Using an analog of results from [32] for CG applied to $AA^{T}y=b$ one can
expect that the size of the term $s_{k}^{T}r_{k+1}$ should be negligible in
comparison to $\|r_{k}\|^{2}$ and that the identity

$\|x-x_{k}\|^{2}-\|x-x_{k+1}\|^{2}=\gamma_{k}\|r_{k}\|^{2}$

holds numerically, until the level of maximal attainable accuracy is reached.
Therefore, for CGNE we obtain an analog to (2) in the form,

$\|x-x_{{\ell}}\|^{2}=\sum_{j={\ell}}^{{k}}\gamma_{j}\|r_{j}\|^{2}+\|x-x_{{k+1}}\|^{2}$

giving the error estimator

$\Delta^{\mathrm{CGNE}}_{{\ell:k}}\equiv\sum_{j={\ell}}^{{k}}\gamma_{j}\|r_{j}\|^{2}\approx\|x-x_{{\ell}}\|^{2}.$

###  3.3 Error estimate in CRAIG

Considering the formula (21) for the iteration index $k-1$, and using
$x_{k}=A^{T}y_{k}$ we obtain

$\left\|x-x_{k-1}\right\|^{2}-\left\|x-x_{k}\right\|^{2}=\|x_{k}-x_{k-1}\|^{2}+2\left(y_{k}-y_{k-1}\right)^{T}A\left(x-x_{k}\right).$

Let us first realize that

$\|x_{k}-x_{k-1}\|^{2}=\zeta_{k}^{2}\|v_{k}\|^{2}=\zeta_{k}^{2}.$

Using a similar technique as in Lemma 1 one can prove that

$A(x-x_{k})=AA^{T}(y-y_{k})=-\zeta_{k}\beta_{k+1}u_{k+1},$

so that, from line 10 of Algorithm 6,

$\left(y_{k}-y_{k-1}\right)^{T}A\left(x-x_{k}\right)=-\zeta_{k}^{2}\beta_{k+1}d_{k}^{T}u_{k+1}.$

Therefore, assuming that the maximal level of accuracy has not been reached
yet and that local orthogonality between $d_{k}$ (direction vector in CG for
(19)) and $u_{k+1}$ (scaled residual vector in CG for (19)) is well preserved,
the identity

$\|x-x_{k-1}\|^{2}-\|x-x_{k}\|^{2}=\zeta_{k}^{2}$

leading to

$\|x-x_{{\ell}}\|^{2}=\sum_{j={\ell}}^{{k}}\zeta^{2}_{j+1}+\|x-x_{{k+1}}\|^{2}$

holds numerically. A detailed rounding error analysis concerning the
preservation of local orthogonality in CRAIG is beyond the scope of this
paper. Finally, the error estimator for CRAIG has the form

$\Delta^{\mathrm{CRAIG}}_{{\ell:k}}\equiv\sum_{j={\ell}}^{{k}}\zeta^{2}_{j+1}\approx\|x-x_{{\ell}}\|^{2}.$

##  4 Error estimation in preconditioned algorithms

In this section we present the preconditioned variants of the algorithms and
derive the associated error estimates. We restrict ourselves to the class of
split preconditioners for normal equations, i.e., the preconditioners formally
transforming the systems (4) and (19) to

$\big{(}L^{-1}A^{T}AL^{-T}\big{)}\big{(}L^{T}x\big{)}=L^{-1}A^{T}b\quad\mathrm{and}\quad\big{(}L^{-1}AA^{T}L^{-T}\big{)}\big{(}L^{T}y\big{)}=L^{-1}b,\
\mbox{respectively}.$

In [6, Section 1], several preconditioning techniques for least-squares
problems are discussed. Most of them can be represented as a split
preconditioner for normal equations. In general, typical representatives of
split preconditioners are incomplete factorizations of matrices $A^{T}A$, or
$AA^{T}$. Efficient codes for computing such factorizations without explicitly
forming $A^{T}A$, or $AA^{T}$ are available; see, e.g., HSL_MI35 from [1].

First we will discuss preconditioning for CGLS and LSQR, and then for CGNE and
CRAIG.

###  4.1 Preconditioned CGLS and LSQR

Let a nonsingular matrix $L\in\mathbb{R}^{n\times n}$ be given. For least-
squares problems, we consider the modification of the original problem (3) in
the form

$x=\operatorname*{arg\,min}_{z\in\mathbb{R}^{n}}\|b-\underbrace{AL^{-T}}_{\hat{A}}\underbrace{L^{T}z}_{\phantom{\hat{I}}\hat{z}\phantom{\hat{I}}}\|\;=\;\operatorname*{arg\,min}_{\hat{z}\in\mathbb{R}^{n}}\|b-\hat{A}\hat{z}\|,$
(22)

leading to the corresponding system of normal equations

$\underbrace{L^{-1}A^{T}}_{{\hat{A}}^{T}}\underbrace{AL^{-T}}_{\hat{A}}\underbrace{L^{T}x}_{\phantom{\hat{I}}\hat{x}\phantom{\hat{I}}}=\underbrace{L^{-1}A^{T}}_{{\hat{A}}^{T}}b.$
(23)

Hence, $L$ can be seen as a split preconditioner for the matrix $A^{T}A$. Let
$\hat{x}_{k}$ be the $k$th CG approximate solution for the preconditioned
system, and define by $x_{k}=L^{-T}\hat{x}_{k}$ an approximation to the
solution $x$ of the original problem (3). Recalling from (23) that
$\hat{x}=L^{T}x$ and $\hat{A}=AL^{-T}$, we have

$\|\hat{x}-\hat{x}_{k}\|^{2}_{\hat{A}^{T}\hat{A}}=(\hat{x}-\hat{x}_{k})^{T}\hat{A}^{T}\hat{A}(\hat{x}-\hat{x}_{k})=(x-x_{k})^{T}LL^{-1}A^{T}AL^{-T}L^{T}(x-x_{k})=\|{x}-{x}_{k}\|^{2}_{{A}^{T}{A}}.$
(24)

The CGLS algorithm for solving (22), based on Algorithm 3, is given in
Algorithm 7. Simple algebraic manipulations show that the residual
$r_{k}=b-Ax_{k}$ for $x_{k}=L^{-T}\hat{x}_{k}$ satisfies $r_{k}=\hat{r}_{k}$.
This allows us to write a variant of preconditioned CGLS as in Algorithm 8,
where the approximations $x_{k}$ and the associated residuals $r_{k}$ are
explicitly computed.

Analogously to Section 2, we can show that the errors satisfy

$\|\hat{x}-\hat{x}_{k}\|^{2}_{\hat{A}^{T}\hat{A}}=\sum_{j={\ell}}^{{k}}\hat{\gamma}_{j}\|\hat{s}_{j}\|^{2}+\|\hat{x}-\hat{x}_{{k+1}}\|^{2}_{\hat{A}^{T}\hat{A}}.$

Therefore we can estimate the error in preconditioned CGLS using the estimator

$\Delta^{\mathrm{PCGLS}}_{{\ell:k}}\equiv\sum_{j={\ell}}^{{k}}\hat{\gamma}_{j}\|\hat{s}_{j}\|^{2}\approx\|x-x_{{\ell}}\|^{2}_{A^{T}A},$

with the same favourable properties as $\Delta^{\mathrm{CGLS}}_{{\ell:k}}$.

The algorithm and the estimator can easily be modified also for the case when
a preconditioner is not available in the factorized form. The _factorization-
free_ variant of PCGLS can be found, e.g., in [29, Algorithm 4].

Algorithm 7 Formally preconditioned CGLS

1: input $\hat{A}$, $\hat{b}$, $\hat{x}_{0}$

2: $\hat{r}_{0}=\hat{b}-\hat{A}\hat{x}_{0}$

3: $\hat{s}_{0}=\hat{p}_{0}=\hat{A}^{T}\hat{r}_{0}$

4: for $k=0,1,\dots$ do

5: $\hat{q}_{k}=\hat{A}\hat{p}_{k}$

6: $\hat{\gamma}_{k}=\|\hat{s}_{k}\|^{2}/\|\hat{q}_{k}\|^{2}$

7: $\hat{x}_{k+1}=\hat{x}_{k}+\hat{\gamma}_{k}\hat{p}_{k}$

8: $\hat{r}_{k+1}=\hat{r}_{k}-\hat{\gamma}_{k}\hat{q}_{k}$

9: $\hat{s}_{k+1}=\hat{A}^{T}\hat{r}_{k+1}$

10: $\hat{\delta}_{k+1}=\|\hat{s}_{k+1}\|^{2}/\|\hat{s}_{k}\|^{2}$

11: $\hat{p}_{k+1}=\hat{s}_{k+1}+\hat{\delta}_{k+1}\hat{p}_{k}$

12: end for

Algorithm 8 Preconditioned CGLS

1: input $A$, $b$, $x_{0}$, $L$

2: $r_{0}={b}-{A}{x}_{0}$

3: $\hat{s}_{0}=\hat{p}_{0}=L^{-1}{A}^{T}{r}_{0}$

4: for $k=0,1,\dots$ do

5: $\hat{t}_{k}=L^{-T}\hat{p}_{k}$

6: $\hat{q}_{k}={A}t_{k}$

7: $\hat{\gamma}_{k}=\|\hat{s}_{k}\|^{2}/\|\hat{q}_{k}\|^{2}$

8: ${x}_{k+1}={x}_{k}+\hat{\gamma}_{k}\hat{t}_{k}$

9: ${r}_{k+1}={r}_{k}-\hat{\gamma}_{k}\hat{q}_{k}$

10: $\hat{s}_{k+1}=L^{-1}{A}^{T}{r}_{k+1}$

11: $\hat{\delta}_{k+1}=\|\hat{s}_{k+1}\|^{2}/\|\hat{s}_{k}\|^{2}$

12: $\hat{p}_{k+1}=\hat{s}_{k+1}+\hat{\delta}_{k+1}\hat{p}_{k}$

13: end for

For LSQR we can proceed analogously to CGLS and apply LSQR directly to (22).
After stopping iterations with $\hat{x}_{k}$ in hands, the approximation to
the solution $x$ of the original system is computed as
$x_{k}=L^{-T}\hat{x}_{k}$. The resulting algorithm is given in Algorithm 9.

From the definition of the estimator $\Delta^{\mathrm{LSQR}}_{{\ell:k}}$
(considered for (22)) and (24), we derive the estimator for preconditioned
LSQR

$\Delta^{\mathrm{PLSQR}}_{{\ell:k}}\equiv\sum_{j={\ell}}^{{k}}\hat{\phi}_{j+1}^{\,2}\approx\|\hat{x}-\hat{x}_{{\ell}}\|^{2}_{\hat{A}^{T}\hat{A}}=\|{x}-{x}_{{\ell}}\|^{2}_{{A}^{T}{A}}.$

Algorithm 9 Preconditioned LSQR

1: input $A$, $b$, $L$

2: $\hat{\beta}_{1}\hat{u}_{1}=b$ $\hat{\beta}_{i}$ are the normalization
coefficients to have $\|\hat{u}_{i}\|=1$, $i=1,2,\ldots$

3: $\hat{\alpha}_{1}\hat{v}_{1}=L^{-1}A^{T}\hat{u}_{1}$ $\hat{\alpha}_{i}$ are
the normalization coefficients to have $\|\hat{v}_{i}\|=1$, $i=1,2,\ldots$

4: $\hat{w}_{1}=\hat{v}_{1}$

5: $\hat{x}_{0}$ = 0

6: $\hat{\bar{\phi}}_{1}=\hat{\beta}_{1}$

7: $\hat{\bar{\rho}}_{1}=\hat{\alpha}_{1}$

8: for $k=1,2,\dots$ do

9:
$\hat{\beta}_{k+1}\hat{u}_{k+1}=\hat{A}\hat{v}_{k}-\hat{\alpha}_{k}\hat{u}_{k}=AL^{-T}\hat{v}_{k}-\hat{\alpha}_{k}\hat{u}_{k}$

10:
$\hat{\alpha}_{k+1}\hat{v}_{k+1}=\hat{A}^{T}\hat{u}_{k+1}-\hat{\beta}_{k+1}\hat{v}_{k}=L^{-1}{A}^{T}\hat{u}_{k+1}-\hat{\beta}_{k+1}\hat{v}_{k}$

11: $\hat{\rho}_{k}=(\hat{\bar{\rho}}_{k}^{2}+\hat{\beta}^{2}_{k+1})^{1/2}$

12: $\hat{c}_{k}=\hat{\bar{\rho}}_{k}/\hat{\rho}_{k}$

13: $\hat{s}_{k}=\hat{\beta}_{k+1}/\hat{\rho}_{k}$

14: $\hat{\theta}_{k+1}=\hat{s}_{k}\hat{\alpha}_{k+1}$

15: $\hat{\bar{\rho}}_{k+1}=-\hat{c}_{k}\hat{\alpha}_{k+1}$

16: $\hat{\phi}_{k}=\hat{c}_{k}\hat{\bar{\phi}}_{k}$

17: $\hat{\bar{\phi}}_{k+1}=\hat{s}_{k}\hat{\bar{\phi}}_{k}$

18: $\hat{x}_{k}=\hat{x}_{k-1}+(\hat{\phi}_{k}/\hat{\rho}_{k})\hat{w}_{k}$

19:
$\hat{w}_{k+1}=\hat{v}_{k+1}+(\hat{\theta}_{k+1}/\hat{\rho}_{k})\hat{w}_{k}$

20: end for

21: $x_{k}=L^{-T}\hat{x}_{k}$

Factorization-free variant of preconditioned LSQR can be found, e.g., in [4,
Algorithm 2].

###  4.2 Preconditioned CGNE and CRAIG

Let a nonsingular $L\in\mathbb{R}^{m\times m}$ be given. We consider the
following modification of the original problem (18)

$x=\operatorname*{arg\,min}_{z\in\mathbb{R}^{n}}\|z\|\quad\mbox{subject
to}\quad\underbrace{L^{-1}A}_{\hat{A}}z=\underbrace{L^{-1}b}_{\hat{b}}\,,$

and, as before, the approximate solutions $x_{k}$ can be obtained when
applying the considered algorithms to the underlying preconditioned system

$\underbrace{L^{-1}A}_{\hat{A}}\underbrace{A^{T}L^{-T}}_{\hat{A}^{T}}\underbrace{L^{T}y}_{\hat{y}}=\underbrace{L^{-1}b}_{\hat{b}}.$

Therefore, $L$ can be seen as a split preconditioner for the matrix $AA^{T}$.
We note that no transformation of the computed approximation is needed to get
an approximation to the solution $x$. Analogously to Section 3.1, we can show
that the preconditioned variants minimize the Euclidean norm of the error
$\|x-x_{k}\|^{2}$, now over the affine space
$x_{0}+\mathcal{K}_{k}(L^{-1}A^{T}AL^{-T},A^{T}L^{-T}L^{-1}r_{0})$.

Algorithm 10 Preconditioned CGNE

1: input ${A}$, ${b}$, ${x}_{0}$, $L$

2: $\hat{r}_{0}=L^{-1}({b}-{A}x_{0})$

3: $\hat{p}_{0}={A}^{T}L^{-T}\hat{r}_{0}$

4: for $k=0,1,2,\dots$ do

5: $\hat{\gamma}{}_{k}={\|\hat{r}_{k}\|^{2}}/{\|\hat{p}_{k}\|^{2}}$

6: $x_{k+1}=x_{k}+\hat{\gamma}{}_{k}\hat{p}_{k}$

7: $\hat{r}_{k+1}=\hat{r}_{k}-\hat{\gamma}{}_{k}L^{-1}{A}\hat{p}_{k}$

8: $\hat{\delta}_{k+1}={\|\hat{r}_{k+1}\|^{2}}/{\|\hat{r}_{k}\|^{2}}$

9: $\hat{p}_{k+1}={A}^{T}L^{-T}\hat{r}_{k+1}+\hat{\delta}_{k+1}\hat{p}_{k}$

10: end for

For CGNE (Algorithm 5) we obtain the preconditioned variant, Algorithm 10.
Whenever the residuals $r_{k}=b-Ax_{k}$ are needed within the iterations, we
can use the transformation $\hat{r}_{k}=L^{-1}r_{k}$ and replace lines 2 and 7
in Algorithm 10 by

$\displaystyle 2:{r}_{0}={b}-{A}x_{0};\quad\hat{r}_{0}=L^{-1}r_{0}$
$\displaystyle
7:{r}_{k+1}={r}_{k}-\hat{\gamma}{}_{k}{A}\hat{p}_{k};\quad\hat{r}_{k}=L^{-1}r_{k}$

We can similarly precondition the GRAIG algorithm (Algorithm 6). A version
without computing the $y_{k}$ vectors is given in Algorithm 11.

Algorithm 11 Preconditioned CRAIG

1: input $A$, $b$, $L$

2: $\hat{\zeta}_{0}=-1$

3: $x_{0}=0$

4: $\hat{\beta}_{1}\hat{u}_{1}=L^{-1}b$ $\hat{\beta}_{i}$ are the
normalization coefficients to have $\|\hat{u}_{i}\|=1$, $i=1,2,\ldots$

5: $\hat{\alpha}_{1}v_{1}=A^{T}L^{-T}\hat{u}_{1}$ $\hat{\alpha}_{i}$ are the
normalization coefficients to have $\|\hat{v}_{i}\|=1$, $i=1,2,\ldots$

6: for $k=1,2,\dots$ do

7:
$\hat{\zeta}_{k}=-\frac{\hat{\beta}_{k}}{\hat{\alpha}_{k}}\hat{\zeta}_{k-1}$

8: $x_{k}=x_{k-1}+\hat{\zeta}_{k}\hat{v}_{k}$

9:
$\hat{\beta}_{k+1}\hat{u}_{k+1}=L^{-1}A\hat{v}_{k}-\hat{\alpha}_{k}\hat{u}_{k}$

10:
$\hat{\alpha}_{k+1}\hat{v}_{k+1}=A^{T}L^{-T}\hat{u}_{k+1}-\hat{\beta}_{k+1}\hat{v}_{k}$

11: end for

Proceeding as in Section 3.2 and Section 3.3, we obtain the error estimators
in PCGNE and PCRAIG,

$\Delta^{\mathrm{PCGNE}}_{{\ell:k}}\equiv\sum_{j={\ell}}^{{k}}\hat{\gamma}_{j}\|\hat{r}_{j}\|^{2}\approx\|x-x_{{\ell}}\|^{2},\quad\Delta^{\mathrm{PCRAIG}}_{{\ell:k}}\equiv\sum_{j={\ell}}^{{k}}\hat{\zeta}^{2}_{j+1}\approx\|x-x_{{\ell}}\|^{2}.$

##  5 Numerical experiments

For the numerical tests we consider several matrices and systems from the
SuiteSparse111https://sparse.tamu.edu, [10] matrix collection, with the sizes

illc1033*: | 1033 | $\times$ | 320
---|---|---|---
illc1850*: | 1850 | $\times$ | 712
well1033*: | 1033 | $\times$ | 320
well1850*: | 1850 | $\times$ | 712
sls: | 1 748 122 | $\times$ | 62 729
Delor64K: | 64 719 | $\times$ | 1 785 345
---|---|---|---
Delor338K: | 343 236 | $\times$ | 887 058
flower_7_4: | 27 693 | $\times$ | 67 593
cat_ears_4_4: | 19 020 | $\times$ | 44 448
lp_pilot: | 1441 | $\times$ | 4860

The matrices marked with an asterisk come together with the right-hand side
$b$. For the remaining problems, we generated $b$ as follows (in MATLAB
notation)

        x = ones(size(A,2),1);
        x(2:2:end) = -2;
        x(5:5:end) = 0;
        b_LN = A*x;
        b_LS = b_LN + randn(size(b_LN))*norm(b_LN);

where $b_{\mathrm{LN}}$ is used for least-norm problems and $b_{\mathrm{LS}}$
for least-squares. The exact solution is then computed using MATLAB build-in
function `lsqminnorm`.

As in [24], we plot the error quantity together with the (adaptive) lower
bound and compare the adaptive value of delay $k-\ell$ with the _ideal_ value,
i.e. the minimal delay ensuring the prescribed accuracy of the bound.

###  5.1 Least-squares problems

In Figures 2–6 we plot the results for our adaptive error estimate in least-
squares problems solved by CGLS and LSQR. We can see that the estimate mostly
follows the error very tightly with nearly the optimal delay, even in the
cases with almost stagnation where the optimal delay is very large; cf. Figure
2. We observe some underestimation in initial iterations for illc1033 matrix
but in later iterations, where one typically needs an estimate for stopping
the solver, the delay is close to the ideal value. For the matrix sls, we note
that the adaptively chosen value of delay is higher than needed.

matrix illc1033, CGLS

iterations

matrix illc1033, LSQR

iterations

Figure 2: Matrix illc1033, CGLS (left) and LSQR (right): error
$\|x-x_{k}\|_{A^{T}A}$ and adaptive error estimate (top), adaptively chosen
delay $k-\ell$ and its ideal value (bottom)

matrix illc1850, CGLS

iterations

matrix illc1850, LSQR

iterations

Figure 3: Matrix illc1850, CGLS (left) and LSQR (right): error
$\|x-x_{k}\|_{A^{T}A}$ and adaptive error estimate (top), adaptively chosen
delay $k-\ell$ and its ideal value (bottom)

matrix well1033, CGLS

iterations

matrix well1033, LSQR

iterations

Figure 4: Matrix well1033, CGLS (left) and LSQR (right): error
$\|x-x_{k}\|_{A^{T}A}$ and adaptive error estimate (top), adaptively chosen
delay $k-\ell$ and its ideal value (bottom)

matrix well1850, CGLS

iterations

matrix well1850, LSQR

iterations

Figure 5: Matrix well1850, CGLS (left) and LSQR (right): error
$\|x-x_{k}\|_{A^{T}A}$ and adaptive error estimate (top), adaptively chosen
delay $k-\ell$ and its ideal value (bottom)

matrix sls, CGLS

iterations

matrix sls, LSQR

iterations

Figure 6: Matrix sls, CGLS (left) and LSQR (right): error
$\|x-x_{k}\|_{A^{T}A}$ and adaptive error estimate (top), adaptively chosen
delay $k-\ell$ and its ideal value (bottom)

###  5.2 Least-norm problems

The results for least-norm problems solved by CGNE and CRAIG are given in
Figures 7–11. Here we observe very satisfactory behaviour with slightly higher
delay for Delor338K and lp_pilot; nevertheless, the adaptively chosen delay
nicely follows the increases and decreases of the ideal value.

matrix Delor64K, CGNE

iterations

matrix Delor64K, CRAIG

iterations

Figure 7: Matrix Delor64K, CGNE (left) and CRAIG (right): error $\|x-x_{k}\|$
and adaptive error estimate (top), adaptively chosen delay $k-\ell$ and its
ideal value (bottom)

matrix Delor338K, CGNE

iterations

matrix Delor338K, CRAIG

iterations

Figure 8: Matrix Delor338K, CGNE (left) and CRAIG (right): error $\|x-x_{k}\|$
and adaptive error estimate (top), adaptively chosen delay $k-\ell$ and its
ideal value (bottom)

matrix flower_7_4, CGNE

iterations

matrix flower_7_4, CRAIG

iterations

Figure 9: Matrix flower_7_4, CGNE (left) and CRAIG (right): error
$\|x-x_{k}\|$ and adaptive error estimate (top), adaptively chosen delay
$k-\ell$ and its ideal value (bottom)

matrix cat_ears_4_4, CGNE

iterations

matrix cat_ears_4_4, CRAIG

iterations

Figure 10: Matrix cat_ears_4_4, CGNE (left) and CRAIG (right): error
$\|x-x_{k}\|$ and adaptive error estimate (top), adaptively chosen delay
$k-\ell$ and its ideal value (bottom)

matrix lp_pilot, CGNE

iterations

matrix lp_pilot, CRAIG

iterations

Figure 11: Matrix lp_pilot, CGNE (left) and CRAIG (right): error $\|x-x_{k}\|$
and adaptive error estimate (top), adaptively chosen delay $k-\ell$ and its
ideal value (bottom)

###  5.3 Preconditioned least-squares problems

In Figures 12–16 we plot the results for our adaptive error estimate in
preconditioned least-squares problems solved by CGLS and LSQR. The
preconditioner $L$ is constructed from $A$ using MATLAB interface of HSL_MI35
from [1], i.e., computing the incomplete Cholesky decomposition of $A^{T}A$,
without explicitly forming it. In the figures, we observe a very satisfactory
behavior of the adaptive error estimate.

matrix illc1033, PCGLS

iterations

matrix illc1033, PLSQR

iterations

Figure 12: Matrix illc1033 with preconditioning, PCGLS (left) and PLSQR
(right): error $\|x-x_{k}\|_{A^{T}A}$ and adaptive error estimate (top),
adaptively chosen delay $k-\ell$ and its ideal value (bottom)

matrix illc1850, PCGLS

iterations

matrix illc1850, PLSQR

iterations

Figure 13: Matrix illc1850 with preconditioning, PCGLS (left) and PLSQR
(right): error $\|x-x_{k}\|_{A^{T}A}$ and adaptive error estimate (top),
adaptively chosen delay $k-\ell$ and its ideal value (bottom)

matrix well1033, PCGLS

iterations

matrix well1033, PLSQR

iterations

Figure 14: Matrix well1033 with preconditioning, PCGLS (left) and PLSQR
(right): error $\|x-x_{k}\|_{A^{T}A}$ and adaptive error estimate (top),
adaptively chosen delay $k-\ell$ and its ideal value (bottom)

matrix well1850, PCGLS

iterations

matrix well1850, PLSQR

iterations

Figure 15: Matrix well1850 with preconditioning, PCGLS (left) and PLSQR
(right): error $\|x-x_{k}\|_{A^{T}A}$ and adaptive error estimate (top),
adaptively chosen delay $k-\ell$ and its ideal value (bottom)

matrix sls, PCGLS

iterations

matrix sls, PLSQR

iterations

Figure 16: Matrix sls with preconditioning, PCGLS (left) and PLSQR (right):
error $\|x-x_{k}\|_{A^{T}A}$ and adaptive error estimate (top), adaptively
chosen delay $k-\ell$ and its ideal value (bottom)

###  5.4 Preconditioned least-norm problems

The results for preconditioned least-norm problems solved by CGNE and CRAIG
are given in Figures 17–21. As in the previous subsection, HSL_MI35 is used to
construct the preconditioner $L$. The experiments confirm that the estimate
can be reliably used also for preconditioned least-norm problems.

matrix Delor64K, PCGNE

iterations

matrix Delor64K, PCRAIG

iterations

Figure 17: Matrix Delor64K with preconditioning, PCGNE (left) and PCRAIG
(right): error $\|x-x_{k}\|$ and adaptive error estimate (top), adaptively
chosen delay $k-\ell$ and its ideal value (bottom)

matrix Delor338K, PCGNE

iterations

matrix Delor338K, PCRAIG

iterations

Figure 18: Matrix Delor338K with preconditioning, PCGNE (left) and PCRAIG
(right): error $\|x-x_{k}\|$ and adaptive error estimate (top), adaptively
chosen delay $k-\ell$ and its ideal value (bottom)

matrix flower_7_4, PCGNE

iterations

matrix flower_7_4, PCRAIG

iterations

Figure 19: Matrix flower_7_4 with preconditioning, PCGNE (left) and PCRAIG
(right): error $\|x-x_{k}\|$ and adaptive error estimate (top), adaptively
chosen delay $k-\ell$ and its ideal value (bottom)

matrix cat_ears_4_4, PCGNE

iterations

matrix cat_ears_4_4, PCRAIG

iterations

Figure 20: Matrix cat_ears_4_4 with preconditioning, PCGNE (left) and PCRAIG
(right): error $\|x-x_{k}\|$ and adaptive error estimate (top), adaptively
chosen delay $k-\ell$ and its ideal value (bottom)

matrix lp_pilot, PCGNE

iterations

matrix lp_pilot, PCRAIG

iterations

Figure 21: Matrix lp_pilot with preconditioning, PCGNE (left) and PCRAIG
(right): error $\|x-x_{k}\|$ and adaptive error estimate (top), adaptively
chosen delay $k-\ell$ and its ideal value (bottom)

##  6 Concluding discussion

In this paper we focused on accurate estimation of errors in CG-like
algorithms for solving least-squares and least-norm problems. Such estimates
are computed just from the available coefficients, and their evaluation is
very cheap. Until a level of maximal attainable accuracy is reached, the
estimates are numerically reliable under the assumption that the local
orthogonality among consecutive vectors is preserved. The consecutive vectors
we have in mind correspond to the underlying CG direction vector with the
iteration index $k-1$ and residual vector with index $k$. In this paper we did
not analyze in detail the preservation of local orthogonality in the
individual algorithms. However, based on previous results [32, 33, 23] for CG
one can expect that such an analysis is doable for all the algorithms
discussed in this paper with analogous results. If needed, for example to
verify the reliability of estimates in a particular application, the local-
orthogonality terms can be evaluated, typically for the price of one extra
inner product.

Our aim in this paper was to obtain an error estimate with a prescribed
relative accuracy $\tau$. At the current iteration $k$, we estimate the
quantity of interest related to some previous iteration $\ell\leq k$. We
developed a heuristic strategy based on our previous results on CG [24] to
make the delay $k-\ell$ as short as possible. Our numerical results show that
the suggested heuristic strategy is robust and reliable, and that the proposed
delay $k-\ell$ is often almost optimal.

The suggested approach provides estimates that represent lower bounds on the
quantity of interest. Nevertheless, once a lower bound $\Delta_{\ell:k}$ with
a prescribed relative accuracy $\tau$ is obtained, the quantity

$\frac{\Delta_{\ell:k}}{1-\tau}$

represents an upper bound; see the discussion in [24, Section 3.1]. Hence, we
can also easily obtain tight (but not guaranteed) upper bounds.

In summary, in CGLS and LSQR we obtain tight estimates of the quantity

$\|x-x_{k}\|_{A^{T}A}^{2}=\left\|r_{k}\right\|^{2}-\left\|b-b_{|\mathcal{R}(A)}\right\|^{2}$

that can be used in stopping criteria of the algorithms. An example of such a
stopping criterion is discussed in [7] and [19]. It is suggested to stop the
iterations when

$\|x-x_{k}\|_{A^{T}A}^{2}\leq\alpha\|A\|\|x_{k}\|+\beta\|b\|,$

where $0\leq\alpha,\beta\ll 1$ are some prescribed tolerances. In CGNE and
CRAIG we are able to estimate efficiently the quantity $\|x-x_{k}\|^{2}.$
Finally, assuming $x_{0}=0,$ our techniques can be straightforwardly applied
for estimating the relative quantities

$\frac{\|x-x_{k}\|_{A^{T}A}^{2}}{\|x\|_{A^{T}A}^{2}}\quad\mbox{and}\quad\frac{\|x-x_{k}\|^{2}}{\|x\|^{2}},$

in least-squares and least-norm problems, respectively; see also [33].

We hope that the results presented in this paper will prove to be useful in
practical computations. They allow to approximate the errors at a negligible
cost during iterations of the considered algorithm, while taking into account
the prescribed relative accuracy of the estimates. The MATLAB codes of the
algorithms with error estimates are available from the GitHub repository [28].

#### Acknowledgement:

The work of Jan Papež has been supported by the Czech Academy of Sciences (RVO
67985840) and by the Grant Agency of the Czech Republic (grant no. 23-06159S).
The authors would like to thank Gérard Meurant for careful reading of the
manuscript and helpful comments, which have greatly improved the presentation.

## References

  * [1] HSL. A collection of Fortran codes for large scale scientific computation. http://www.hsl.rl.ac.uk/.
  * [2] M. Arioli, Generalized Golub–Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., 34 (2013), pp. 571–592.
  * [3] M. Arioli and S. Gratton, Linear regression models, least-squares problems, normal equations, and stopping criteria for the conjugate gradient method, Comput. Phys. Commun., 183 (2012), pp. 2322–2336.
  * [4] S. R. Arridge, M. M. Betcke, and L. Harhanen, Iterated preconditioned LSQR method for inverse problems on unstructured grids, Inverse Problems, 30 (2014), pp. 075009, 27.
  * [5] Å. Björck, T. Elfving, and Z. Strakoš, Stability of conjugate gradient and Lanczos methods for linear least squares problems, SIAM J. Matrix Anal. Appl., 19 (1998), pp. 720–736.
  * [6] R. Bru, J. Marín, J. Mas, and M. Tůma, Preconditioned iterative methods for solving linear least squares problems, SIAM J. Sci. Comput., 36 (2014), pp. A2002–A2022.
  * [7] X.-W. Chang, C. C. Paige, and D. Titley-Peloquin, Stopping criteria for the iterative solution of linear least squares problems, SIAM J. Matrix Anal. Appl., 31 (2009), pp. 831–852.
  * [8] E. J. Craig, Iteration procedures for simultaneous equations, PhD thesis, Massachusetts Institute of Technology, 1954.
  * [9] , The $N$-step iteration procedures, J. Math. and Phys., 34 (1955), pp. 64–73.
  * [10] T. A. Davis and Y. Hu, The university of Florida sparse matrix collection, ACM Trans. Math. Softw., 38 (2011).
  * [11] R. Estrin, D. Orban, and M. A. Saunders, Euclidean-norm error bounds for SYMMLQ and CG, SIAM J. Matrix Anal. Appl., 40 (2019), pp. 235–253.
  * [12] , LNLQ: an iterative method for least-norm problems with an error minimization property, SIAM J. Matrix Anal. Appl., 40 (2019), pp. 1102–1124.
  * [13] , LSLQ: an iterative method for linear least-squares with an error minimization property, SIAM J. Matrix Anal. Appl., 40 (2019), pp. 254–275.
  * [14] D. K. Faddeev and V. N. Faddeeva, Computational Methods of Linear Algebra, W. H. Freeman and Co., San Francisco, 1963.
  * [15] G. H. Golub and Z. Strakoš, Estimates in quadratic formulas, Numer. Algorithms, 8 (1994), pp. 241–268.
  * [16] A. Greenbaum, Behavior of slightly perturbed Lanczos and conjugate-gradient recurrences, Linear Algebra Appl., 113 (1989), pp. 7–63.
  * [17] E. Hallman, Sharp 2-norm error bounds for LSQR and the conjugate gradient method, SIAM J. Matrix Anal. Appl., 41 (2020), pp. 1183–1207.
  * [18] M. R. Hestenes and E. Stiefel, Methods of conjugate gradients for solving linear systems, J. Research Nat. Bur. Standards, 49 (1952), pp. 409–436.
  * [19] P. Jiránek and D. Titley-Peloquin, Estimating the backward error in LSQR, SIAM J. Matrix Anal. Appl., 31 (2009/10), pp. 2055–2074.
  * [20] E. Kaasschieter, Preconditioned conjugate gradients for solving singular systems, J. Comput. Appl. Math., 24 (1988), pp. 265–275.
  * [21] J. Liesen and Z. Strakoš, Krylov Subspace Methods, Numerical Mathematics and Scientific Computation, Oxford University Press, Oxford, 2013\.
  * [22] G. Meurant, Estimates of the $l_{2}$ norm of the error in the conjugate gradient algorithm, Numer. Algorithms, 40 (2005), pp. 157–169.
  * [23] G. Meurant, The Lanczos and Conjugate Gradient Algorithms, from Theory to Finite Precision Computations, SIAM, Philadelphia, 2006.
  * [24] G. Meurant, J. Papež, and P. Tichý, Accurate error estimation in CG, Numer. Algorithms, 88 (2021), pp. 1337–1359.
  * [25] G. Meurant and P. Tichý, On computing quadrature-based bounds for the A-norm of the error in conjugate gradients, Numer. Algorithms, 62 (2013), pp. 163–191.
  * [26] C. C. Paige, Bidiagonalization of matrices and solutions of the linear equations, SIAM J. Numer. Anal., 11 (1974), pp. 197–209.
  * [27] C. C. Paige and M. A. Saunders, LSQR: an algorithm for sparse linear equations and sparse least squares, ACM Trans. Math. Software, 8 (1982), pp. 43–71.
  * [28] J. Papež and P. Tichý, CG-like methods with error estimate, GitHub repository.  https://github.com/JanPapez/CGlike-methods-with-error-estimate, 2023\.
  * [29] S. Regev and M. A. Saunders, SSAI: A symmetric sparse approximate inverse preconditioner for the conjugate gradient methods PCG and PCGLS, 2020\.  Unpublished report, available at https://web.stanford.edu/group/SOL/reports/20SSAI.pdf.
  * [30] M. A. Saunders, Solution of sparse rectangular systems using LSQR and CRAIG, BIT, 35 (1995), pp. 588–604.
  * [31] G. L. G. Sleijpen, H. A. van der Vorst, and D. R. Fokkema, ${\rm BiCGstab}(l)$ and other hybrid Bi-CG methods, Numer. Algorithms, 7 (1994), pp. 75–109.
  * [32] Z. Strakoš and P. Tichý, On error estimation in the conjugate gradient method and why it works in finite precision computations, Electron. Trans. Numer. Anal., 13 (2002), pp. 56–80.
  * [33] , Error estimation in preconditioned conjugate gradients, BIT, 45 (2005), pp. 789–817.