# A comparison of mixed precision iterative refinement approaches for least-
squares problems

Erin Carson Faculty of Mathematics and Physics, Charles University, {carson,
<EMAIL_ADDRESS>Both authors were supported by Charles
University Research Centre program No. UNCE/24/SCI/005, the Exascale Computing
Project (17-SC-20-SC), a collaborative effort of the U.S. Department of Energy
Office of Science and the National Nuclear Security Administration, and by the
European Union (ERC, inEXASCALE, 101075632). Views and opinions expressed are
those of the authors only and do not necessarily reflect those of the European
Union or the European Research Council. Neither the European Union nor the
granting authority can be held responsible for them. Ieva
Daužickaitė11footnotemark: 1

###### Abstract

Various approaches to iterative refinement (IR) for least-squares problems
have been proposed in the literature and it may not be clear which approach is
suitable for a given problem. We consider three approaches to IR for least-
squares problems when two precisions are used and review their theoretical
guarantees, known shortcomings and when the method can be expected to
recognize that the correct solution has been found, and extend uniform
precision analysis for an IR approach based on the semi-normal equations to
the two-precision case. We focus on the situation where it is desired to
refine the solution to the working precision level. It is shown that the IR
methods exhibit different sensitivities to the conditioning of the problem and
the size of the least-squares residual, which should be taken into account
when choosing the IR approach. We also discuss a new approach that is based on
solving multiple least-squares problems.

###### keywords:

mixed precision, iterative refinement, least-squares, semi-normal equations

###### AMS:

65F05, 65F20

##  1 Introduction

Iterative refinement is used to reduce the impact of finite precision errors
in computations. When solving linear systems of equations $Ax=b$, the solution
can be refined by repeatedly computing the residual in a possibly higher
precision and solving the system with the residual on the right-hand side to
obtain an update to the solution, that is, computing $x_{0}$ in some way and
repeating the following steps:

  1. 1.

Compute $r_{i}=b-Ax_{i}$.

  2. 2.

Solve $A\delta x_{i}=r_{i}$.

  3. 3.

Update $x_{i+1}=x_{i}+\delta x_{i}$.

Various approaches for solving the correction equation in step 2 with
convergence guarantees depending on the condition number of the coefficient
matrix $A$ can be used, e.g. [2, 10].

More care may be needed in the case of least-squares problems

$\min_{x}\|b-Ax\|_{2},$ (1)

where $A\in\mathbb{R}^{m\times n}$, $m\geq n$, is full rank,
$x\in\mathbb{R}^{n}$, and $b\in\mathbb{R}^{m}$. The least-squares residual
$r=b-Ax$ may be non-zero and this invalidates some of the convergence
guarantees that are available in the linear system case.

There is a renewed interest in iterative refinement for least-squares problems
[14], and laterally in employing low precision computations for such problems
[12, 19]. In this note, we thus examine three approaches to least-squares
iterative refinement (LSIR) in two precisions proposed in the literature, show
how two of the approaches can be combined so that the solution updates come
from solving a few least-squares problems independently, and discuss their
(lack of) convergence guarantees. We also investigate when the methods can be
guaranteed to recognize that a solution correct to a required accuracy is
found in the sense that the new computed update does not corrupt the solution.
Further, the forward error analysis for the semi-normal equations LSIR
approach is extended from the uniform precision case to the case where two
precisions are used. In this case, we determine the limiting accuracy and the
convergence conditions. In the original form, LSIR strategies employ a direct
solve for computing the update and are based on a QR decomposition of $A$ or
its $R$ factor. This may be appropriate when $A$ is dense, but iterative
methods may be preferred for sparse $A$. We include some numerical experiments
to illustrate the tendencies when iterative methods are used, but theoretical
convergence guarantees are not provided.

The discussed LSIR procedures were originally introduced assuming the data are
in single precision and the computations (except computing the residual) are
performed in single precision as well [4, 16]. The higher precision for the
residual computations is achieved by accumulating the inner products in double
precision. We present general versions of these LSIR procedures in two
precisions with unit roundoffs $u$ (working precision) and $u_{r}$ (residual
computations) with $u_{r}\leq u^{2}$.

The focus of our convergence analysis is the relative forward error and when
it can be guaranteed to reach the working precision level $\mathcal{O}(u)$.
One could argue that solving ill-conditioned problems to high accuracy is an
unnecessary task; note, however, that if the working precision is set to a low
precision then the notion of an ill-conditioned problem changes accordingly.
For example, when $u$ is set to IEEE half precision then problems with
condition number as low as $10^{3}$ can be considered to be ill-conditioned in
the working precision even though it is a moderately conditioned problem in
the usual double precision. See Table 1 for parameters of some IEEE precision
formats.

Arithmetic | $x_{min}$ | $x_{max}$ | $u$
---|---|---|---
fp16 (half) | $6.1\times 10^{-5}$ | $6.6\times 10^{4}$ | $4.9\times 10^{-4}$
fp32 (single) | $1.2\times 10^{-38}$ | $3.4\times 10^{38}$ | $6.0\times 10^{-8}$
fp64 (double) | $2.2\times 10^{-308}$ | $1.8\times 10^{308}$ | $1.1\times 10^{-16}$
Table 1: Parameters for IEEE floating point arithmetics: the smallest positive
number $x_{min}$, the largest positive number $x_{max}$, and unit roundoff
$u$.

This note is structured as follows. We discuss the least-squares system
approach to LSIR in Section 2. The semi-normal equations and augmented system
approaches follow in Sections 3 and 4, respectively. A short discussion on
using iterative solvers for LSIR is presented in Section 5. We comment on how
to combine the least-squares system and augmented system approaches in Section
6. Numerical experiments illustrate the theory in Section 7 and we summarize
and give insight into choosing a solver in practice in Section 8.

###  1.1 Notation

In the following, hats are used to refer to the computed quantities; for
example, $\widehat{r}$ is the computed version of $r$. We use a standard model
of floating point arithmetic and define

$\tilde{\gamma}_{n}=\frac{cnu}{1-cnu}\quad\text{and}\quad\tilde{\gamma}_{n}^{(r)}=\frac{cnu_{r}}{1-cnu_{r}},$

where $c$ is a small constant independent of $n$. When $c=1$ we use
$\gamma_{n}$ or $\gamma_{n}^{(r)}$. We use $\lesssim$ when we omit higher
order terms which are insignificant in the expression. $\kappa(A)$ denotes the
$L_{2}$ norm condition number of $A$, that is,
$\kappa(A)=\|A\|_{2}\|A^{\dagger}\|_{2}$, where $A^{\dagger}$ is the Moore-
Penrose pseudoinverse of $A$.

##  2 Least-squares system approach

Golub [15] proposed a straightforward extension of the IR procedure to least-
squares problems, which we specify in two precisions in Algorithm 2.1. We call
it the LS approach.

Algorithm 2.1 Two-precision LSIR based on the least-squares system approach

Input: $A\in\mathbb{R}^{m\times n}$, $b\in\mathbb{R}^{m}$, initial solution
$x_{0}\in\mathbb{R}^{n}$, precisions $u_{r}$ and $u$, where $u_{r}\leq u^{2}$
Output: approximate solution $x$

1:$i=0$

2:while not converged do

3: Compute $r_{i}=b-Ax_{i}$ $\triangleright$ $u_{r}$

4: Solve $\min_{\delta x_{i}}\|r_{i}-A\delta x_{i}\|_{2}$ via QR decomposition
$\triangleright$ $u$

5: Update $x_{i+1}=x_{i}+\delta x_{i}$ $\triangleright$ $u$

6: $i=i+1$

7:end while

No convergence guarantees are provided by Golub, but it is noted that this
LSIR procedure is likely to converge only if the initial approximation $x_{0}$
is of sufficient accuracy. Businger and Golub [9] implement this procedure and
point out that “there is no assurance, however, that all digits of the final
solution will be correct.” The theoretical analysis by Golub and Wilkinson
[16] reveals some shortcomings of the approach in Algorithm 2.1. It is shown
that when Householder QR decomposition is used to solve the correction
equation in step 4, the correct solution may not be recognized (as described
below) unless the system is nearly compatible, that is, $\|r\|_{2}$ is nearly
zero. In the following, we restate and explain their argument in the general
case, i.e., assuming that the solve in step 4 is performed in some precision
$u$. The error in computing the residual $r_{i}$ is ignored.

First, we note that the correction $\delta\widehat{x}_{i}$ computed in step 4
is the exact solution of a perturbed problem (see, e.g., [17, Theorem 20.3])

$\displaystyle\min_{\delta x_{i}}\|r_{i}+\Delta r_{i}-(A+\Delta
A_{i})\delta\widehat{x}_{i}\|_{2},\textrm{ where}$ $\displaystyle\|\Delta
A_{i}\|_{F}\leq\tilde{\gamma}_{mn}\|A\|_{F}\textrm{ and }\|\Delta
r_{i}\|_{2}\leq\tilde{\gamma}_{mn}\|r_{i}\|_{2}.$

Then using [17, Theorem 20.1] under the assumption that
$\tilde{\gamma}_{mn}\kappa(A)<1$, the error between the true correction
$\delta x_{i}$ and the computed $\delta\widehat{x}_{i}$ is bounded as

$\|\delta\widehat{x}_{i}-\delta
x_{i}\|_{2}\leq\frac{\tilde{\gamma}_{mn}}{1-\tilde{\gamma}_{mn}\kappa(A)}\left(\kappa(A)\|\delta
x_{i}\|_{2}+\kappa(A)^{2}\frac{\|\delta r_{i}\|_{2}}{\|A\|_{2}}\right),$ (2)

where $\delta r_{i}=r_{i}-A\delta x_{i}$.

To illustrate the possible failing of Algorithm 2.1, assume that
$x_{i}=x^{*}$, where $x^{*}$ is the minimizer of (1). In this case, we have
$r_{i}=b-Ax^{*}=r^{*}$, and thus $\delta x_{i}=0$ and $\delta r_{i}=r^{*}$.
Ideally in this case, we would have
$\|\delta\widehat{x}_{i}\|_{2}/\|x^{*}\|_{2}=\mathcal{O}(u)$ for the
_computed_ $\delta\widehat{x}_{i}$, since we would like to recognize that we
have at least found the solution to working accuracy and thus the refinement
should be terminated at this point. Unfortunately, from (2), we can only say
that

$\|\delta\widehat{x}_{i}\|_{2}\leq\frac{\tilde{\gamma}_{mn}}{1-\tilde{\gamma}_{mn}\kappa(A)}\kappa(A)^{2}\frac{\|r^{*}\|_{2}}{\|A\|_{2}}.$

This means that in order to guarantee that
$\|\delta\widehat{x}_{i}\|_{2}/\|x^{*}\|_{2}=O(u)$ is achievable, we must have
that

$\kappa(A)^{2}\frac{\|r^{*}\|_{2}}{\|A\|_{2}\|x^{*}\|_{2}}<1.$ (3)

When the least-squares problem is almost compatible, that is, $\|r^{*}\|_{2}$
is small, the condition (3) is satisfied for a wider range of coefficient
matrices $A$. Note that the condition is independent of $u$. Golub and
Wilkinson [16] thus conclude that “whatever precision of computation is used
there will be right-hand sides for which iterative refinement will never give
solutions which are correct to working accuracy.”

##  3 Semi-normal equations approach

Golub and Wilkinson [16] briefly discuss an approach to iterative refinement
that employs the semi-normal equations; the approach is attributed to Kahan.
We present it in Algorithm 3.1.

Algorithm 3.1 Two-precision LSIR based on the semi-normal equations approach

Input: $A\in\mathbb{R}^{m\times n}$, $b\in\mathbb{R}^{m}$,
$R\in\mathbb{R}^{n\times n}$ factor of the Householder QR decomposition of
$A$, initial solution $x_{0}\in\mathbb{R}^{n}$, precisions $u_{r}$ and $u$,
where $u_{r}\leq u^{2}$
Output: approximate solution $x$

1:$i=0$

2:while not converged do

3: Compute $r_{i}=b-Ax_{i}$ $\triangleright$ $u_{r}$

4: Solve $R^{T}R\delta x_{i}=A^{T}r_{i}$ $\triangleright$ Computing
$A^{T}r_{i}$ in $u_{r}$, solving in $u$

5: Update $x_{i+1}=x_{i}+\delta x_{i}$ $\triangleright$ $u$

6: $i=i+1$

7:end while

The update is computed by solving two triangular systems via forward and
backward substitution. Golub and Wilkinson explore when this approach
recognizes the correct solution $x^{*}$ _if_ it is found. They consider the
case when $r_{i}$ and $A^{T}r_{i}$ are computed in higher (double) precision
$u_{r}$ and only the computed $A^{T}r_{i}$ is cast to the lower (single)
precision $u$. Note that if $x_{i}=x^{*}$ and the error in computing
$fl(A^{T}r_{i})$ is ignored, then the right-hand side of the correction
equation in step 4 of Algorithm 3.1 is zero and thus $\delta x_{i}=0$. If
$x_{i}=x^{*}$, $\|A\|_{2}=1$ and we account for the error in computing
$fl(A^{T}r_{i})=A^{T}r^{*}+A^{T}h_{i}=A^{T}h_{i}$ in $u_{r}$, it is shown that
the true solution $\delta\widetilde{x}_{i}$ to the correction equation

$\widehat{R}^{T}\widehat{R}\delta\widetilde{x}_{i}=A^{T}h_{i}$

satisfies

$\|\delta\widetilde{x}_{i}\|_{2}\leq u_{r}\kappa(A)^{2}\|r^{*}\|_{2}.$

We can write this in a form similar to (3). If we set $u_{r}=u^{2}$, then we
can guarantee $\|\delta\widetilde{x}_{i}\|_{2}/\|x^{*}\|_{2}=\mathcal{O}(u)$
when

$\kappa(A)^{2}\frac{\|r^{*}\|_{2}}{\|x^{*}\|_{2}}<u^{-1}.$ (4)

A bound for the computed solution $\delta\widehat{x}_{i}$ (where the errors in
the triangular solves are also accounted for) is not provided in [16]. The
forward error $\|\delta\widetilde{x}_{i}-\delta\widehat{x}_{i}\|_{2}$ for two
triangular solves with $\widehat{R}$ and $\widehat{R}^{T}$ in precision $u$
can be bounded using

$\delta\widehat{x}_{i}=(\widehat{R}+\Delta R_{2})^{-1}(\widehat{R}^{T}+\Delta
R_{1})^{-1}A^{T}h_{i},$

where $\|\Delta R_{i}\|_{2}\leq n^{1/2}\gamma_{n}\|\widehat{R}\|_{2}$, and
$A\widehat{R}^{-1}=\tilde{Q}-E\widehat{R}^{-1}$, where $\widehat{R}$ is
computed via a Householder QR algorithm, $\tilde{Q}^{T}\tilde{Q}=I$ and
$\|E\|_{2}\leq n^{1/2}\widetilde{\gamma}_{mn}\|A\|_{2}$; see, for example,
[17, Theorems 8.5 and 19.4]. Then employing $\|h\|_{2}\leq u_{r}\|r^{*}\|_{2}$
and ignoring higher order terms gives

$\|\delta\widetilde{x}_{i}-\delta\widehat{x}_{i}\|_{2}\lesssim
n^{1/2}\gamma_{n}u_{r}\kappa(A)^{2}\|r^{*}\|_{2}.$

We can guarantee
$\|\delta\widetilde{x}_{i}-\delta\widehat{x}_{i}\|_{\infty}/\|x^{*}\|_{2}=\mathcal{O}(u)$
if

$\kappa(A)^{2}\frac{\|r^{*}\|_{2}}{\|x^{*}\|_{2}}<u_{r}^{-1},$

which clearly holds whenever (4) holds. Note that unlike (3) the condition (4)
depends on the precision used and thus problems caused by large residual or
ill-conditioning can be addressed by increasing the precision.

In [6], Björck derives a bound on the forward error for what he calls the
“corrected” semi-normal equations approach, in which one step of iterative
refinement is performed and a uniform precision is used throughout the
computation. His goal is to obtain a solution with forward error equivalent to
the error given by a backward stable LS solver rather than obtain a forward
error to the level $\mathcal{O}(u)$. Recently Epperly [14] extended Björck’s
analysis to the case when $R$ is computed via a QR factorization $\Omega
A=QR$, where $\Omega$ is a random sketching matrix, and $x_{0}$ is obtained as
$x_{0}=R^{-1}Q^{T}\Omega b$; if a high quality $R$ is obtained, then the
method is forward stable. In the following, we analyze the convergence of the
forward error in the semi-normal equations approach to LSIR, adapting the
analysis for the corrected semi-normal equations in uniform precision of [6]
to the mixed precision case described in Algorithm 3.1.

###  3.1 The forward error

We first consider the computed factor $\widehat{R}$. It is assumed that
$\widehat{R}$ is computed via the backward stable Householder QR algorithm in
precision with unit roundoff $u$ and hence we can write (see, e.g., [17,
Theorem 19.4])

$\displaystyle A+E=\tilde{Q}\widehat{R},\textrm{ where }$ (5)
$\displaystyle\|E\|_{F}\leq\tilde{\gamma}_{mn}\|A\|_{F},$ (6)

and $\tilde{Q}$ is orthogonal. It follows that

$\|A\|_{2}=\|\tilde{Q}\widehat{R}-E\|_{2}\leq\|\tilde{Q}\widehat{R}\|_{2}+\|E\|_{2}\leq\|\widehat{R}\|_{2}+\tilde{\gamma}_{mn}\|A\|_{F}\leq\|\widehat{R}\|_{2}+n^{1/2}\tilde{\gamma}_{mn}\|A\|_{2},$

and assuming that $n^{1/2}\tilde{\gamma}_{mn}<1$, we thus have

$\|A\|_{2}\leq\underbrace{\left(1-n^{1/2}\tilde{\gamma}_{mn}\right)^{-1}}_{\alpha}\|\widehat{R}\|_{2}.$
(7)

We also assume that $\kappa(A)$ and $u$ are such that

$c_{1}(m,n)u\kappa(A)<1,$

where $c_{1}(m,n)$ is a constant depending on $m$ and $n$. Then $\widehat{R}$
is invertible and

$\kappa(\widehat{R})=\kappa(A)\left(1+\mathcal{O}(u)\kappa(A)\right).$ (8)

Our goal in this section is to bound the error
$\|x^{*}-\widehat{x}_{i+1}\|_{2}$. We do this by accounting for the finite
precision error in steps 3 to 5 of Algorithm 3.1 as follows:

$\displaystyle\widehat{r}_{i}=b-A\widehat{x}_{i}+f,$ (9)
$\displaystyle\widehat{R}^{T}\widehat{R}\delta\widehat{x}_{i}=A^{T}\widehat{r}_{i}+k,$
(10) $\displaystyle\widehat{x}_{i+1}=\widehat{x}_{i}+\delta\widehat{x}_{i}+e,$
(11)

where the vector $f$ accounts for the finite precision error in computing
$fl(r_{i})$, $k$ accounts for the error in computing
$fl(A^{T}\widehat{r}_{i})$ and the triangular solves with $\widehat{R}$, and
$e$ accounts for the error in adding $\widehat{x}_{i}$ and
$\delta\widehat{x}_{i}$. Using (9)-(11), we write

$\displaystyle x^{*}-\widehat{x}_{i+1}=$
$\displaystyle\,x^{*}-\widehat{x}_{i}-\delta\widehat{x}_{i}-e$
$\displaystyle=$
$\displaystyle\,x^{*}-\widehat{x}_{i}-(\widehat{R}^{T}\widehat{R})^{-1}\left(A^{T}\widehat{r}_{i}+k\right)-e$
$\displaystyle=$
$\displaystyle\,x^{*}-\widehat{x}_{i}-(\widehat{R}^{T}\widehat{R})^{-1}\left(A^{T}(b-A\widehat{x}_{i}+f)+k\right)-e$
$\displaystyle=$
$\displaystyle\,M(x^{*}-\widehat{x}_{i})-(\widehat{R}^{T}\widehat{R})^{-1}\left(A^{T}f+k\right)-e,$
(12)

where $M=I-(\widehat{R}^{T}\widehat{R})^{-1}A^{T}A$ and we used
$A^{T}b=A^{T}Ax$ in the last equation. We first aim to bound $\|M\|_{2}$ and
for this we define

$\bar{Q}=A\widehat{R}^{-1}.$

Using (5) - (8), we have

$\displaystyle\bar{Q}=\tilde{Q}-E\widehat{R}^{-1},$
$\displaystyle\|\bar{Q}\|_{2}\leq\|\tilde{Q}\|_{2}+\|E\widehat{R}^{-1}\|_{2}\leq
1+n^{1/2}\tilde{\gamma}_{mn}\|A\|_{2}\|\widehat{R}^{-1}\|_{2}\lesssim
1+n^{1/2}\alpha\tilde{\gamma}_{mn}\kappa(A).$

We then write

$M=\widehat{R}^{-1}\left(I-\bar{Q}^{T}\bar{Q}\right)\widehat{R}.$ (13)

Note that

$I-\bar{Q}^{T}\bar{Q}=\tilde{Q}^{T}E\widehat{R}^{-1}+\widehat{R}^{-T}E^{T}\tilde{Q}-\widehat{R}^{-T}E^{T}E\widehat{R}^{-1}.$

Taking the norm and using (6) and (7) gives

$\displaystyle\|I-\bar{Q}^{T}\bar{Q}\|_{2}\leq$
$\displaystyle\,2\|E\|_{2}\|\widehat{R}^{-1}\|_{2}+\|E\|_{2}^{2}\|\widehat{R}^{-1}\|_{2}^{2}$
$\displaystyle\leq$
$\displaystyle\,2n^{1/2}\tilde{\gamma}_{mn}\|A\|_{2}\|\widehat{R}^{-1}\|_{2}+n\tilde{\gamma}_{mn}^{2}\|A\|_{2}^{2}\|\widehat{R}^{-1}\|_{2}^{2}$
$\displaystyle\leq$
$\displaystyle\,2n^{1/2}\tilde{\gamma}_{mn}\alpha\kappa(\widehat{R})+n\tilde{\gamma}_{mn}^{2}\alpha^{2}\kappa(\widehat{R})^{2}.$

Combining this with (13) and (8) results in

$\displaystyle\|M\|_{2}\leq$
$\displaystyle\,2n^{1/2}\tilde{\gamma}_{mn}\alpha\kappa(\widehat{R})^{2}+n\tilde{\gamma}_{mn}^{2}\alpha^{2}\kappa(\widehat{R})^{3}$
$\displaystyle\leq$
$\displaystyle\,2n^{1/2}\tilde{\gamma}_{mn}\alpha\kappa(A)^{2}\left(1+\mathcal{O}(u\kappa(A))\right)^{2}+n\tilde{\gamma}_{mn}^{2}\alpha^{2}\kappa(A)^{3}\left(1+\mathcal{O}(u\kappa(A))\right)^{3}$
$\displaystyle\lesssim$
$\displaystyle\,2n^{1/2}\tilde{\gamma}_{mn}\alpha\kappa(A)^{2}.$ (14)

We proceed by bounding the norm of
$(\widehat{R}^{T}\widehat{R})^{-1}\left(A^{T}f+k\right)$. The error $f$ can be
bounded as

$\|f\|_{2}\leq
n^{1/2}\tilde{\gamma}_{n}^{(r)}\left(\|A\|_{2}\|\widehat{x}_{i}\|_{2}+\|b\|_{2}\right).$
(15)

We write the error vector $k$ as $k=k_{1}+k_{2}$, where $k_{1}$ accounts for
the error in computing $fl(A^{T}\widehat{r}_{i})$ and $k_{2}$ accounts for the
error in two triangular solves. Then

$\|k_{1}\|_{2}\leq n^{1/2}\gamma_{n}^{(r)}\|A\|_{2}\|\widehat{r}\|_{2}\lesssim
n^{1/2}\gamma_{n}^{(r)}\|A\|_{2}\left(\|A\|_{2}\|\widehat{x}_{i}\|_{2}+\|b\|_{2}\right),$
(16)

where we ignore the terms in $u_{r}^{2}$. To bound $\|k_{2}\|_{2}$, we use
standard results on the backward error in triangular solves, e.g. [17, Theorem
8.5], and obtain

$\left(\widehat{R}^{T}+\Delta R_{1}^{T}\right)\left(\widehat{R}+\Delta
R_{2}\right)\delta\widehat{x}_{i}=A^{T}\widehat{r}_{i}+k_{1},$ (17)

where $\|\Delta R_{j}\|_{F}\leq\gamma_{n}\|\widehat{R}\|_{F}$ for $j=1,2$.
Then using (10) we have

$k+\left(\widehat{R}^{T}\Delta R_{2}+\Delta R_{1}^{T}\widehat{R}+\Delta
R_{2}\Delta R_{1}^{T}\right)\delta\widehat{x}_{i}=k_{1}$

and

$k_{2}=k-k_{1}=-\left(\widehat{R}^{T}\Delta R_{2}+\Delta
R_{1}^{T}\widehat{R}+\Delta R_{2}\Delta
R_{1}^{T}\right)\delta\widehat{x}_{i}.$

Multiplying with $\widehat{R}^{-T}$ on the left and moving $\widehat{R}$ gives

$\widehat{R}^{-T}k_{2}=-\left(\Delta
R_{2}\widehat{R}^{-1}+\widehat{R}^{-T}\Delta R_{1}^{T}+\widehat{R}^{-T}\Delta
R_{2}\Delta R_{1}^{T}\widehat{R}^{-1}\right)\widehat{R}\delta\widehat{x}_{i}.$

We thus have

$\|\widehat{R}^{-T}k_{2}\|_{2}\leq\left(2n^{1/2}\gamma_{n}\kappa(\widehat{R})+n\gamma_{n}^{2}\kappa(\widehat{R})^{2}\right)\|\widehat{R}\delta\widehat{x}_{i}\|_{2}$
(18)

and need to bound $\|\widehat{R}\delta\widehat{x}_{i}\|_{2}$. To do this, we
write (17) as

$\widehat{R}^{T}\left(I+\widehat{R}^{-T}\Delta R_{1}^{T}\right)\left(I+\Delta
R_{2}\widehat{R}^{-1}\right)\widehat{R}\delta\widehat{x}_{i}=A^{T}\widehat{r}_{i}+k_{1}.$

Assuming that $\|\widehat{R}^{-T}\Delta R_{1}^{T}\|_{2}<1$ (which is
essentially equivalent to $\kappa(A)<u^{-1}$), we use the approximation
$\left(I+\widehat{R}^{-T}\Delta R_{1}^{T}\right)^{-1}\approx
I-\widehat{R}^{-T}\Delta R_{1}^{T}$ and an equivalent argument for
$\left(I+\Delta R_{2}\widehat{R}^{-1}\right)^{-1}$ to obtain

$\widehat{R}\delta\widehat{x}_{i}\approx\left(I-\Delta
R_{2}\widehat{R}^{-1}\right)\left(I-\widehat{R}^{-T}\Delta
R_{1}^{T}\right)\left(\widehat{R}^{-T}A^{T}\widehat{r}_{i}+\widehat{R}^{-T}k_{1}\right).$

Taking the norm and using $A^{T}r^{*}=0$ gives

$\displaystyle\|\widehat{R}\delta\widehat{x}_{i}\|_{2}\leq$
$\displaystyle\,\left(1+n^{1/2}\gamma_{n}\kappa(\widehat{R})\right)^{2}\left(\|\widehat{R}^{-T}A^{T}\widehat{r}_{i}\|_{2}+\|\widehat{R}^{-T}k_{1}\|_{2}\right)$
$\displaystyle=$
$\displaystyle\,\left(1+n^{1/2}\gamma_{n}\kappa(\widehat{R})\right)^{2}\left(\|\bar{Q}^{T}(\widehat{r}_{i}-r^{*})\|_{2}+\|\widehat{R}^{-T}k_{1}\|_{2}\right).$

Combining this with (18), (16) and ignoring terms in
$u^{2}\kappa(\widehat{R})^{2}$ and $uu_{r}\kappa(\widehat{R})^{2}$ gives

$\displaystyle\|\widehat{R}^{-T}k_{2}\|_{2}\lesssim$
$\displaystyle\,2n^{1/2}\gamma_{n}\kappa(\widehat{R})\|\bar{Q}^{T}\|_{2}\|\widehat{r}_{i}-r^{*}\|_{2}$
$\displaystyle\lesssim$
$\displaystyle\,2n^{1/2}\gamma_{n}\kappa(\widehat{R})\|r^{*}-(b-A\widehat{x}_{i})\|_{2}.$
(19)

Now using (12), (14), (15), (16), (19) and $\|e\|_{2}\leq
u\|\widehat{x}_{i+1}\|_{2}$ [17, eq. 2.5] we obtain

$\displaystyle\|x^{*}-\widehat{x}_{i+1}\|_{2}\leq$
$\displaystyle\,\|M(x^{*}-\widehat{x}_{i})\|_{2}+\|(\widehat{R}^{T}\widehat{R})^{-1}\left(A^{T}f+k\right)\|_{2}+\|e\|_{2}$
$\displaystyle\leq$
$\displaystyle\,\|M\|_{2}\|x^{*}-\widehat{x}_{i}\|_{2}+\|\widehat{R}^{-1}\bar{Q}^{T}f\|_{2}+\|\widehat{R}^{-1}\widehat{R}^{-T}k\|_{2}+\|e\|_{2}$
$\displaystyle\lesssim$
$\displaystyle\,2n^{1/2}\tilde{\gamma}_{mn}\alpha\kappa(A)^{2}\|x^{*}-\widehat{x}_{i}\|_{2}$
$\displaystyle\,+n^{1/2}\tilde{\gamma}_{n}^{(r)}\|\widehat{R}^{-1}\|_{2}\left(\|A\|_{2}\|\widehat{x}_{i}\|_{2}+\|b\|_{2}\right)$
$\displaystyle\,+n^{1/2}\gamma_{n}^{(r)}\|\widehat{R}^{-1}\|_{2}^{2}\|A\|_{2}\left(\|A\|_{2}\|\widehat{x}_{i}\|_{2}+\|b\|_{2}\right)$
$\displaystyle\,+2n^{1/2}\gamma_{n}\kappa(\widehat{R})\|\widehat{R}^{-1}\|_{2}\|r^{*}-(b-A\widehat{x}_{i})\|_{2}$
$\displaystyle\,+u\|\widehat{x}_{i+1}\|_{2}$ $\displaystyle\leq$
$\displaystyle\,2n^{1/2}\tilde{\gamma}_{mn}\alpha\kappa(A)^{2}\|x^{*}-\widehat{x}_{i}\|_{2}$
$\displaystyle\,+2n^{1/2}\gamma_{n}\alpha\kappa(A)^{2}\frac{\|r^{*}-(b-A\widehat{x}_{i})\|_{2}}{\|A\|_{2}}$
$\displaystyle\,+u\|\widehat{x}_{i+1}\|_{2}+n^{1/2}\gamma_{n}^{(r)}\kappa(A)^{2}\frac{\|A\|_{2}\|\widehat{x}_{i}\|_{2}+\|b\|_{2}}{\|A\|_{2}}.$

Using $r^{*}-(b-A\widehat{x}_{i})=A(\widehat{x}_{i}-x^{*})$ gives

$\displaystyle\|x^{*}-\widehat{x}_{i+1}\|_{2}\leq$
$\displaystyle\,n^{1/2}\tilde{\gamma}_{mn}\alpha\kappa(A)^{2}\|x^{*}-\widehat{x}_{i}\|_{2}$
$\displaystyle\,+u\|\widehat{x}_{i+1}\|_{2}+n^{1/2}\gamma_{n}^{(r)}\kappa(A)^{2}\frac{\|A\|_{2}\|\widehat{x}_{i}\|_{2}+\|b\|_{2}}{\|A\|_{2}}.$
(20)

We can thus conclude that the solution is improved with every LSIR iteration
if $n^{1/2}\tilde{\gamma}_{mn}\alpha\kappa(A)^{2}<1$, which is essentially

$\kappa(A)<u^{-1/2}.$ (21)

Using $\widehat{x}_{i}=\widehat{x}_{i+1}+\mathcal{O}(u)=x^{*}+\mathcal{O}(u)$
in the limit and $\|b\|_{2}\leq\|r^{*}\|_{2}+\|A\|_{2}\|x^{*}\|_{2}$, we can
conclude that the limiting normwise relative error is of the order

$u+n^{1/2}\gamma_{n}^{(r)}\kappa(A)^{2}\left(2+\frac{\|r^{*}\|_{2}}{\|A\|_{2}\|x^{*}\|_{2}}\right).$

We thus achieve $\mathcal{O}(u)$ relative error if

$u_{r}\kappa(A)^{2}\left(2+\frac{\|r^{*}\|_{2}}{\|A\|_{2}\|x^{*}\|_{2}}\right)<u$
(22)

and if $u_{r}=u^{2}$ we require

$\kappa(A)^{2}\left(2+\frac{\|r^{*}\|_{2}}{\|A\|_{2}\|x^{*}\|_{2}}\right)<u^{-1},$

which is satisfied when (21) holds and

$\frac{\|r^{*}\|_{2}}{\|A\|_{2}\|x^{*}\|_{2}}=\mathcal{O}(1).$

Notice that as long as the above holds, $\|r^{*}\|_{2}$ does not affect the
limiting error.

We summarize our results in the following theorem.

######  Theorem 1.

Let $\widehat{x}_{i}$ be the solution computed in iteration $i$ of Algorithm
3.1 and $x^{*}$ be the true solution of (1). Then the relative forward error
after $i$ iterations is bounded as

$\displaystyle\frac{\|x^{*}-\widehat{x}_{i}\|_{2}}{\|x^{*}\|_{2}}\leq$
$\displaystyle\,c_{1}(m,n)u\kappa(A)^{2}\frac{\|x^{*}-\widehat{x}_{i-1}\|_{2}}{\|x^{*}\|_{2}}$
$\displaystyle\,+u\frac{\|\widehat{x}_{i}\|_{2}}{\|x^{*}\|_{2}}+c_{2}(m,n)u_{r}\kappa(A)^{2}\frac{\|A\|_{2}\|\widehat{x}_{i-1}\|_{2}+\|b\|_{2}}{\|A\|_{2}\|x^{*}\|_{2}},$

where $c_{1}(m,n)$ and $c_{2}(m,n)$ are constants depending on $m$ and $n$.

We now contrast this with the corrected semi-normal equations method of Björck
[6], which is equivalent to running Algorithm 3.1 with only one refinement
iteration, computing $x_{0}$ by solving the semi-normal equations in precision
$u$, and using $u_{r}=u$. Under these assumptions on $x_{0}$ but keeping
$u_{r}$ general, combining the bound for $\|x^{*}-x_{0}\|_{2}$ in [6, eq.
3.11] with (3.1) gives the following bound for the forward error after one
iteration of LSIR:

$\displaystyle\|x^{*}-\widehat{x}_{1}\|_{2}\leq$
$\displaystyle\,c(m,n)u^{2}\kappa(A)^{3}\left(\|x^{*}\|_{2}+\frac{\|b\|_{2}}{\|A\|_{2}}\right)$
$\displaystyle\,+u\|\widehat{x}_{1}\|_{2}+n^{1/2}\gamma_{n}^{(r)}\kappa(A)^{2}\frac{\|A\|_{2}\|\widehat{x}_{0}\|_{2}+\|b\|_{2}}{\|A\|_{2}}.$

We can thus expect LSIR to reach
$\|x^{*}-\widehat{x}_{1}\|_{2}/\|x^{*}\|_{2}=\mathcal{O}(u)$ and thus converge
after one iteration if $u_{r}=u^{2}$ and

$u\kappa(A)^{3}\left(2+\frac{\|r^{*}\|_{2}}{\|A\|_{2}\|x^{*}\|_{2}}\right)<1,$
(23)

which is equivalent to $\kappa(A)<u^{-1/3}$ when
$\frac{\|r^{*}\|_{2}}{\|A\|_{2}\|x^{*}\|_{2}}=\mathcal{O}(1)$, and
$\kappa(A)\frac{\|r^{*}\|_{2}^{1/3}}{\|A\|_{2}^{1/3}\|x^{*}\|_{2}^{1/3}}<u^{-1/3}$
for large $\|r^{*}\|_{2}$.

Note that the condition (23) is stricter than the condition (22) in the case
that $u_{r}=u^{2}$. If $u_{r}=u$, then the accuracy of Algorithm 3.1 is
limited by the last term in (3.1), which can be written as

$u\kappa(A)^{2}\left(2+\frac{\|r^{*}\|_{2}}{\|A\|_{2}\|x^{*}\|_{2}}\right)$

and we can expect $\|x^{*}-\widehat{x}_{1}\|_{2}/\|x^{*}\|_{2}=\mathcal{O}(u)$
only if $\kappa(A)=\mathcal{O}(1)$ and
$\frac{\|r^{*}\|_{2}}{\|A\|_{2}\|x^{*}\|_{2}}=\mathcal{O}(1)$. In the first
iterations of LSIR however we can expect the first term to dominate and thus
we can get some improvement if (21) holds.

It is thus clear that in many cases, there is a benefit to allowing more than
one iteration of refinement, as mentioned in [7, Chapter 6] for the uniform
precision case. We note that there are no assumptions on how $x_{0}$ is
computed in our general analysis above. If it is of poor quality, LSIR may
indeed need a few iterations to correct the error.

We note that Björck mentions performing fixed precision iterative refinement
with normal equations, where the $R$ factor is obtained via a Cholesky
decomposition of $A^{T}A$ [7, Chapter 6]. Our numerical experiments (not
reported here) show that two precision LSIR with normal equations stops
working for problems with slightly smaller condition number than in
experiments with LSIR with semi-normal equations in Section 7.

##  4 Augmented system approach

Motivated by the sensitivity of the LS system approach to the size of the
residual, Björck [4] exploits the fact that $x$ is the solution of an
overdetermined least-squares problem (1) if and only if it satisfies the
orthogonality property [7, Theorem 1.1.2]

$A^{T}r=0.$

Björck proposes performing iterative refinement on an augmented
$(m+n)\times(m+n)$ linear system of equations

$\underbrace{\begin{bmatrix}I&A\\\
A^{T}&0\end{bmatrix}}_{\widetilde{A}}\begin{bmatrix}r\\\
x\end{bmatrix}=\begin{bmatrix}b\\\ 0\end{bmatrix},$ (24)

which is equivalent to the normal equations. We describe the approach in
Algorithm 4.1; see [13] for a variant of this algorithm and some
implementation details. Note that both $r$ and $x$ are refined. In order to
compare this approach with the LS system approach further, in the following
sections we explore the conditions for recognizing the correct solution
$x_{i}$ and the available convergence guarantees when the correction equation
for the system (24) is solved via a QR factorization. A general stability
analysis when (24) is solved via different QR factorizations is provided in
[8].

Algorithm 4.1 Two-precision LSIR based on the augmented system approach

Input: $A\in\mathbb{R}^{m\times n}$, $b\in\mathbb{R}^{m}$, initial solution
$x_{0}$, precisions $u_{r}$ and $u$ where $u_{r}\leq u^{2}$
Output: approximate solution $x$

1:$i=0$

2:Compute $r_{0}=b-Ax_{0}$ $\triangleright$ $u$

3:while not converged do

4: Compute $\begin{bmatrix}f_{i}\\\ g_{i}\end{bmatrix}=\begin{bmatrix}b\\\
0\end{bmatrix}-\begin{bmatrix}I&A\\\
A^{T}&0\end{bmatrix}\begin{bmatrix}r_{i}\\\ x_{i}\end{bmatrix}$
$\triangleright$ $u_{r}$

5: Solve  $\begin{bmatrix}I&A\\\ A^{T}&0\end{bmatrix}\begin{bmatrix}\delta
r_{i}\\\ \delta x_{i}\end{bmatrix}=\begin{bmatrix}f_{i}\\\ g_{i}\end{bmatrix}$
(25) $\triangleright$ $u$

6: Update $\begin{bmatrix}r_{i+1}\\\
x_{i+1}\end{bmatrix}=\begin{bmatrix}r_{i}\\\
x_{i}\end{bmatrix}+\begin{bmatrix}\delta r_{i}\\\ \delta x_{i}\end{bmatrix}$
$\triangleright$ $u$

7: $i=i+1$

8:end while

We note that Björck discusses other LSIR approaches based on the remainders
$f_{i}$ and $g_{i}$ as in Algorithm 4.1 and using either the semi-normal
equations, where the triangular factor is obtained via QR, or the normal
equations method employing the Cholesky decomposition [5]. A detailed
convergence analysis of these schemes is not provided in [5] and our numerical
experiments, which are not reported here, show that it performs similar to the
semi-normal equations approach described in the previous section. We thus do
not explore this approach further in this note.

###  4.1 Correction to true solution

We wish to understand if the augmented system approach avoids the restrictive
condition (3) for recognizing the correct solution in the LS system approach.
This requires a perturbation result for a system of the form in (24) and
results on the finite precision error incurred when solving such a system via
a QR decomposition. Björck provides the perturbation result, which we restate
here.

######  Theorem 2 ([4], Theorem 1).

Let $\delta\widehat{x}$ and $\delta\widehat{r}$ satisfy a perturbed system of
equations

$\begin{bmatrix}I&A+\Delta A_{1}\\\ A^{T}+\Delta
A_{2}^{T}&0\end{bmatrix}\begin{bmatrix}\delta\widehat{r}\\\
\delta\widehat{x}\end{bmatrix}=\begin{bmatrix}f+\Delta f\\\ g+\Delta
g\end{bmatrix},$

where $A$ is a given $m\times n$ matrix of rank $n$, and let $\delta x$ and
$\delta r$ denote the solution to the corresponding unperturbed system.

Assume that there exists an orthogonal matrix $Q$ such that the matrices
$Q(A+H_{i})$, $i=1,2$ are upper triangular, and let the perturbation satisfy
the bounds

$\displaystyle\|\Delta A_{i}\|_{2}\leq$ $\displaystyle\eta\|A\|_{2},\ i=1,2,$
$\displaystyle\|\Delta f\|_{2}\leq$ $\displaystyle\tau_{2},$
$\displaystyle\|\Delta g\|_{2}\leq$ $\displaystyle\tau_{1}\|A\|_{2}.$

Then, provided

$\beta=(\sqrt{2}+1)\kappa(A)\eta<1,$

it holds that

$\begin{pmatrix}\|\delta\widehat{r}-\delta r\|_{2}\\\
\|A\|_{2}\|\delta\widehat{x}-\delta
x\|_{2}\end{pmatrix}\leq\begin{pmatrix}\chi&1\\\
\chi^{2}&\chi\end{pmatrix}\begin{pmatrix}\eta\|\delta r\|_{2}+\tau_{1}\\\
\eta\|A\|_{2}\|\delta x\|_{2}+\tau_{2}\end{pmatrix},$

where

$\chi=(1-\beta)^{-1/2}\kappa(A).$ (26)

The finite precision error depends on the solution process. Let
$A=Q\begin{pmatrix}R\\\ 0\end{pmatrix}=Q\widetilde{R}$ be the QR decomposition
of $A$, where $Q\in\mathbb{R}^{m\times m}$ and $R\in\mathbb{R}^{n\times n}$.
We can write the block $LDL^{T}$ decomposition of $\widetilde{A}$ as

$\widetilde{A}=\underbrace{\begin{pmatrix}I&0\\\
\widetilde{R}^{T}Q^{T}&I\end{pmatrix}}_{L}\underbrace{\begin{pmatrix}I&0\\\
0&-R^{T}R\end{pmatrix}}_{D}\begin{pmatrix}I&Q\widetilde{R}\\\
0&I\end{pmatrix}.$

Then the solution of (25) is

$\displaystyle\begin{bmatrix}\delta r_{i}\\\ \delta x_{i}\end{bmatrix}=$
$\displaystyle\,L^{-T}D^{-1}L^{-1}\begin{bmatrix}f_{i}\\\ g_{i}\end{bmatrix}$
$\displaystyle=$ $\displaystyle\,\begin{pmatrix}I&-Q\widetilde{R}\\\
0&I\end{pmatrix}\begin{pmatrix}I&0\\\
0&-R^{-1}R^{-T}\end{pmatrix}\begin{pmatrix}I&0\\\
-\widetilde{R}^{T}Q^{T}&I\end{pmatrix}\begin{bmatrix}f_{i}\\\
g_{i}\end{bmatrix}$ $\displaystyle=$
$\displaystyle\,\begin{bmatrix}Q\begin{pmatrix}0&0\\\
0&I\end{pmatrix}Q^{T}f_{i}+Q\begin{pmatrix}R^{-T}\\\ 0\end{pmatrix}g_{i}\\\
\begin{pmatrix}R^{-1}&0\end{pmatrix}Q^{T}f_{i}-R^{-1}R^{-T}g_{i}\end{bmatrix},$
(27)

which can be computed in the following procedure proposed in [4]:

$\displaystyle h=$ $\displaystyle\,R^{-T}g_{i},$ (28) $\displaystyle k=$
$\displaystyle\,Q^{T}f_{i}=\begin{bmatrix}k_{1}\\\ k_{2}\end{bmatrix}$ (29)
$\displaystyle\delta r=$ $\displaystyle\,Q\begin{bmatrix}h\\\
k_{2}\end{bmatrix},$ (30) $\displaystyle\delta x=$
$\displaystyle\,R^{-1}(k_{1}-h).$ (31)

The results of the finite precision analysis are given in the following
theorem.

######  Theorem 3 ([17], Theorem 20.4).

Let $A\in\mathbb{R}^{m\times n}$ be of full rank $n\leq m$ and suppose the
system (25) is solved via a Householder QR factorization of $A$ as described
above in (28)-(31). The computed $\delta\widehat{x}$ and $\delta\widehat{r}$
satisfy

$\begin{bmatrix}I&A+\Delta A_{1}\\\ A^{T}+\Delta
A_{2}^{T}&0\end{bmatrix}\begin{bmatrix}\delta\widehat{r}\\\
\delta\widehat{x}\end{bmatrix}=\begin{bmatrix}f+\Delta f_{1}\\\ g+\Delta
g_{1}\end{bmatrix},$

where

$\displaystyle|\Delta A_{i}|\leq$ $\displaystyle mn\tilde{\gamma}_{m}G|A|,\
i=1,2,$ $\displaystyle|\Delta f_{1}|\leq$ $\displaystyle
m^{1/2}n\tilde{\gamma}_{m}(H_{1}|f|+H_{2}|\delta\widehat{r}|),$
$\displaystyle|\Delta g_{1}|\leq$ $\displaystyle
m^{1/2}n\tilde{\gamma}_{m}|A^{T}|H_{3}|\delta\widehat{r}|,$

where $\|G\|_{F}=\|H_{i}\|_{F}=1$, $i=1:3$.

We can now use Theorems 2 and 3 to obtain bounds on the forward error in both
$\delta\widehat{x}$ and $\delta\widehat{r}$. We set

$\displaystyle\eta=$ $\displaystyle\,mn^{3/2}\tilde{\gamma}_{m},$
$\displaystyle\tau_{2}=$
$\displaystyle\,m^{1/2}n\tilde{\gamma}_{m}\left(\|\widehat{f}\|+\|\delta\widehat{r}\|\right),$
$\displaystyle\tau_{1}=$
$\displaystyle\,m^{1/2}n^{3/2}\tilde{\gamma}_{m}\|\delta\widehat{r}\|_{2}.$

Then if

$(\sqrt{2}+1)mn^{3/2}\tilde{\gamma}_{m}\kappa(A)<1,$

we have

$\displaystyle\|A\|_{2}\|\delta\widehat{x}-\delta x\|_{2}\leq$
$\displaystyle\,\chi^{2}\left(mn^{3/2}\tilde{\gamma}_{m}\|\delta
r\|_{2}+m^{1/2}n^{3/2}\tilde{\gamma}_{m}\|\delta\widehat{r}\|_{2}\right)$
$\displaystyle+\chi\left(mn^{3/2}\tilde{\gamma}_{m}\|A\|_{2}\|\delta
x\|_{2}+m^{1/2}n\tilde{\gamma}_{m}\left(\|\widehat{f}\|_{2}+\|\delta\widehat{r}\|_{2}\right)\right)$
$\displaystyle\leq$
$\displaystyle\,mn^{3/2}\tilde{\gamma}_{m}\chi\left(\chi(\|\delta
r\|_{2}+\|\delta\widehat{r}\|_{2})+\|A\|_{2}\|\delta
x\|_{2}+\|\widehat{f}\|_{2}+\|\delta\widehat{r}\|_{2}\right)$ (32)

and

$\displaystyle\|\delta\widehat{r}-\delta r\|_{2}\leq$
$\displaystyle\,\chi\left(mn^{3/2}\tilde{\gamma}_{m}\|\delta
r\|_{2}+m^{1/2}n^{3/2}\tilde{\gamma}_{m}\|\delta\widehat{r}\|_{2}\right)$
$\displaystyle+\left(mn^{3/2}\tilde{\gamma}_{m}\|A\|_{2}\|\delta
x\|_{2}+m^{1/2}n\tilde{\gamma}_{m}\left(\|\widehat{f}\|_{2}+\|\delta\widehat{r}\|_{2}\right)\right)$
$\displaystyle\leq$
$\displaystyle\,mn^{3/2}\tilde{\gamma}_{m}\left(\chi(\|\delta
r\|_{2}+\|\delta\widehat{r}\|_{2})+\|A\|_{2}\|\delta
x\|_{2}+\|\widehat{f}\|_{2}+\|\delta\widehat{r}\|_{2}\right).$ (33)

We now explore what happens if $x_{i}=x^{*}$. In order to determine what the
ideal update should be, we first ignore the errors in computing $f_{i}$ and
$g_{i}$ and in the solve. Then we have

$\begin{bmatrix}f_{i}\\\ g_{i}\end{bmatrix}=\begin{bmatrix}b\\\
0\end{bmatrix}-\begin{bmatrix}I&A\\\
A^{T}&0\end{bmatrix}\begin{bmatrix}r_{i}\\\
x^{*}_{i}\end{bmatrix}=\begin{bmatrix}r^{*}-r_{i}\\\
-A^{T}r_{i}\end{bmatrix}.$

Combining this with (27) gives

$\begin{bmatrix}\delta r_{i}\\\ \delta
x_{i}\end{bmatrix}=\begin{bmatrix}Q\begin{pmatrix}0&0\\\
0&I\end{pmatrix}Q^{T}(r^{*}-r_{i})-Q\begin{pmatrix}R^{-T}\\\
0\end{pmatrix}A^{T}r_{i}\\\
\begin{pmatrix}R^{-1}&0\end{pmatrix}Q^{T}(r^{*}-r_{i})+R^{-1}R^{-T}A^{T}r_{i}\end{bmatrix}=\begin{bmatrix}r^{*}-r_{i}\\\
0\end{bmatrix},$ (34)

where we use $Q=\begin{bmatrix}Q_{1}&Q_{2}\end{bmatrix}$ with
$Q_{1}\in\mathbb{R}^{m\times n}$ and $Q_{2}\in\mathbb{R}^{m\times(m-n)}$ and
the fact that $Q_{1}^{T}r^{*}=0$. Thus the true solution to the correction
equation is $\delta x_{i}=0$.

We hence want the computed solution to satisfy
$\|\delta\widehat{x}_{i}\|_{2}/\|x^{*}\|_{2}=\mathcal{O}(u)$. Using (32) with
$\delta x_{i}=0$ and the fact that $\chi>1$, the condition for recognizing the
correct solution is

$\frac{\|\delta
r_{i}\|_{2}+2\|\delta\widehat{r}_{i}\|_{2}+\|\widehat{f}_{i}\|_{2}}{\|A\|_{2}\|x^{*}\|_{2}}<\frac{1}{\chi^{2}}.$
(35)

Note that (35) is expected to be significantly less restrictive compared to
(3), because the left-hand side in (35) includes the _corrections_ to the
residual compared to the residual itself in (3).

###  4.2 Initial rate of convergence and sufficient conditions for convergence

Björck analyzes the convergence of LSIR when the augmented system is solved
via the Householder QR approach detailed in the previous section [4, Section
6]. It is assumed that the computations are performed in single precision and
the inner products are accumulated in double precision throughout the
computations except when updating $r_{i+1}$ and $x_{i+1}$. Björck points out
that using higher precision is only essential in the computation of the
residuals and if the accumulation is omitted from other computations, then the
error bounds would increase by at most a factor of $m$.

Björck shows that the initial rate of improvement is proportional to the term

$\rho_{i}=\frac{\|\hat{x}_{i}-x^{*}\|_{2}}{\|\hat{x}_{i-1}-x^{*}\|_{2}}<c(m,n)u\kappa(A),$

where $c(m,n)$ is a constant. Note that the bound does not depend on the size
of the residual and thus it is possible to get some improvement in the first
iterations even for problems with large $\|r^{*}\|_{2}$. The dependency on
$\kappa(A)$ rather than $\kappa(A)^{2}$ also shows that initial improvement is
possible even for highly ill-conditioned problems.

However, $\kappa(A)^{2}$ arises in the conditions that guarantee a small
relative error in the solution. It is shown that

$\lim_{i\rightarrow\infty}\|\widehat{x}_{i}-x^{*}\|_{2}=\mathcal{O}(u)\|x^{*}\|_{2}$

is guaranteed to hold if

$\chi^{2}\frac{\|r^{*}\|_{2}}{\|A\|_{2}\|x^{*}\|_{2}}\leq\frac{1}{{c}_{1}(m,n,u,u_{r})},$
(36)

where ${c}_{1}(m,n,u,u_{r})$ is a constant that depends on $m$, $n$, $u$, and
$u_{r}$. Ignoring the dimensional constants, we can approximate
${c}_{1}(m,n,u,u_{r})$ by its leading term $u$, and also approximate $\chi$ by
$\kappa(A)$ to obtain the simplified convergence condition

$\kappa(A)^{2}\frac{\|r^{*}\|_{2}}{\|A\|_{2}\|x^{*}\|_{2}}\leq u^{-1}.$ (37)

Björck also shows that a sufficient condition for

$\lim_{i\rightarrow\infty}\|\widehat{r}_{i}-r^{*}\|_{2}=\mathcal{O}(u)\|r^{*}\|_{2}$

to hold is

$\frac{\|A\|_{2}\|x^{*}\|_{2}}{\|r^{*}\|_{2}}\leq{c}_{2}(m,n,u,u_{r}),$ (38)

where ${c}_{2}(m,n,u,u_{r})$ is another constant that depends on $m$, $n$,
$u$, and $u_{r}$.

The conditions (37) and (3) are similar, yet in the augmented system approach
it guarantees that LSIR converges in $x$ and the dependency on the unit
roundoff makes it less restrictive, whereas for the least-squares system
approach we can only guarantee that the correct solution will be recognized
_if_ the method finds it, that is, there is no guarantee that the correct
solution will be found.

The following sensitivities of the augmented system approach to LSIR can be
observed.

  * •

The right-hand side of (36) is inversely proportional to $\chi^{2}$ which is
proportional to $\kappa(A)^{2}$. We can thus expect the procedure to become
less reliable in correcting $x$ when $\kappa(A)$ increases.

  * •

The left-hand side of (38) grows when $\|r^{*}\|_{2}$ is reduced. It is thus
suggested that the refinement of the residual is sensitive to the size of it
and if the residual norm is very small then the procedure may fail to correct
it to the required accuracy.

NB: Björck points out that for any matrix $A$ there exist right-hand sides $b$
such that either (36) or (38) does not hold. However, at least one of (36) or
(38) _always_ holds.

###  4.3 Scaling of the augmented system

Theory for the iterative refinement of linear systems of equations can also be
applied to analyze the refinement of (24). When (25) in Algorithm 4.1 is
solved via backward and forward substitution using the $LDL^{T}$ decomposition
of $\widetilde{A}$ computed in precision $u$ or using the QR factors of $A$ as
described in the previous section, then the forward and backward error of the
augmented system is guaranteed to converge if $\kappa(\widetilde{A})$ is
sufficiently less than $u^{-1}$; see, e.g., [17, Chapters 12 and 20.5]. The
convergence is thus independent of the size of the residual, but depends on
$\kappa(\widetilde{A})$ instead of $\kappa(A)$. Björck [4] showed that
$\kappa(\widetilde{A})$ relates to the singular values of $A$ and can be
changed significantly with the scaling

$\begin{pmatrix}I&0\\\ 0&\alpha^{-1}I\end{pmatrix}\begin{pmatrix}I&A\\\
A^{T}&0\end{pmatrix}\begin{pmatrix}\alpha I&0\\\
0&I\end{pmatrix}=\begin{pmatrix}\alpha I&A\\\
A^{T}&0\end{pmatrix}=\widetilde{A}_{\alpha},$

where $\alpha\in\mathbb{R}$ is positive and is not related to $\alpha$ in the
previous section. Then

$\kappa(\widetilde{A}_{\alpha})=\frac{\alpha+\sqrt{\alpha^{2}+4\sigma_{max}(A)^{2}}}{\min\\{2\alpha,\sqrt{\alpha^{2}+4\sigma_{min}(A)^{2}}-\alpha\\}},$

where $\sigma_{max}(A)$ and $\sigma_{min}(A)$ are the largest and smallest
singular values of $A$, respectively. If we set
$\alpha=2^{-1/2}\sigma_{min}(A)$, then $\kappa(\widetilde{A}_{\alpha})$
reaches its minimum value and $\kappa(A)\leq\kappa(\widetilde{A}_{\alpha})\leq
2\kappa(A)$; note that $\sigma_{min}(A)$ is expensive to compute. The case
without scaling, that is, setting $\alpha=1$ gives

$\kappa(\widetilde{A}_{\alpha})=\frac{1+\sqrt{1+4\sigma_{max}(A)^{2}}}{\min\\{2,\sqrt{1+4\sigma_{min}(A)^{2}}-1\\}},$

which can be close to $\kappa(A)^{2}$ if $\sigma_{min}(A)\ll 1$.

Arioli, Duff, and de Rijk [3] note that scaling not only reduces
$\kappa(\widetilde{A})$, but also may improve the quality of the computed
$LDL^{T}$ factorization. They show that for some ill-conditioned problems
scaling with $\alpha$ set to a non-optimal value can enable convergence or
improve the quality of the computed solution when the error does not converge
to the $\mathcal{O}(u)$ level.

##  5 Iterative solvers

The theoretical analysis in the previous sections assumes that the least-
squares problems are solved via Householder QR decomposition. The cost of
computing this decomposition may be prohibitive when $A$ is very large or
sparse. In this case, iterative solvers may be used for the correction
equations.

Irrespective of the solver choice, we are not aware of the existence of
conditions that guarantee convergence for the LS system approach to LSIR, and
as we have already mentioned, Golub and Wilkinson note that the method will
not converge for some right-hand sides. If the LS approach is chosen
regardless, then LSQR, which is a popular iterative method and only accesses
$A$ via a procedure for computing matrix-vector products with $A$ and $A^{T}$,
may be preferred by the user [20].

The semi-normal equations approach is essentially direct and requires a good
quality $R$ factor of the QR decomposition. We note that the analysis can be
applicable when $R$ is computed via methods other than Householder QR as long
as they are backward stable, that is, (5) and (6) hold. Iterative solvers
could, however, be used to find $\delta x_{i}$ via the normal equations
$A^{T}A\delta x_{i}=A^{T}b$ instead of the semi-normal equations in step 4 of
Algorithm 3.1. $A^{T}A$ is a symmetric positive definite matrix and thus the
CG method can be used (see e.g. [21]), although for CG there is no guarantee
of backward stability in the usual sense. Backward stability can be guaranteed
for reasonably well-conditioned $A$ if preconditioned GMRES is used [18].

The augmented system approach with an iterative solver for the correction
equations has been studied before. Various solvers and preconditioners
(possibly computed in a low precision $u_{f}$) can be used as long as the
solution satisfies certain conditions [11]; the authors in [11] use left-
preconditioned GMRES. Scott and Tůma [22] use this approach with split-
preconditioned GMRES and MINRES to solve sparse-dense LS problems. As we
discussed in Section 4.3, convergence theory for iterative refinement for
linear systems of equations can be applied in this case and thus convergence
can be guaranteed depending on the condition number of the (potentially
preconditioned) augmented matrix. The performance of the iterative solver
itself will also depend on the properties of this matrix. We thus note that
computing a good scaling factor and preconditioner is paramount to ensure good
performance of the scheme.

##  6 Combining least-squares and augmented system approaches

We further explore whether both $x$ and $r$ can be updated by solving least-
squares problems, which also allows the use of various iterative solvers. We
do this by considering the inverse of the augmented system matrix

$\widetilde{A}^{-1}=\begin{pmatrix}P_{N}&(A^{\dagger})^{T}\\\
A^{\dagger}&-(A^{T}A)^{-1}\end{pmatrix},$

where $P_{N}=I-AA^{\dagger}$ is an orthogonal projector onto the null-space of
$A^{T}$. The update in (25) is then

$\begin{bmatrix}\delta r_{i}\\\ \delta
x_{i}\end{bmatrix}=\begin{pmatrix}P_{N}f_{i}+(A^{\dagger})^{T}g_{i}\\\
A^{\dagger}f_{i}-(A^{T}A)^{-1}g_{i}\end{pmatrix}.$ (39)

We observe that $(A^{\dagger})^{T}g_{i}$ is a minimum-norm solution to the
under-determined least-squares problem

$\min_{y_{1}}\|g_{i}-A^{T}y_{1}\|_{2}.$ (40)

Note that $A^{\dagger}f_{i}$ is the solution of the least-squares problem

$\min_{y_{2}}\|f_{i}-Ay_{2}\|_{2},$ (41)

and $-(A^{T}A)^{-1}g_{i}=-R^{-1}R^{-T}g_{i}$, where $R$ comes from the QR
decomposition of $A$, is the solution of

$R^{T}Ry_{3}=-g_{i},$

which is equivalent to the semi-normal equation approach to solving

$\min_{y_{3}}\|r_{i}-Ay_{3}\|_{2}$ (42)

when the right-hand side of the semi-normal equations is computed in a higher
precision. Note that, however, since $r_{i}$ is used instead of $A^{T}r_{i}$,
it is likely that the method is sensitive to the size of $\|r\|_{2}$.

We can thus obtain $\delta r_{i}$ and $\delta x_{i}$ in (39) via solving the
least-squares problems (40), (41), and (42), and computing the projection
$P_{N}f_{i}$. Note that the latter can be computed as
$(I-Q_{1}Q_{1}^{T})f_{i}$ or $Q_{2}Q_{2}^{T}f_{i}$, where $Q_{1}$ and $Q_{2}$
are as in (34). We summarize the procedure in Algorithm 6.1. Note that all the
least-squares problems can be solved at the same time and a block least-
squares solver may be employed.

Algorithm 6.1 Two-precision LSIR combining least-squares and augmented system
approaches

Input: $A\in\mathbb{R}^{m\times n}$, $b\in\mathbb{R}^{m}$, initial solution
$x_{0}\in\mathbb{R}^{n}$, projector $P_{N}$ onto the the null-space of
$A^{T}$, precisions $u_{r}$ and $u$ where $u_{r}\leq u$
Output: approximate solution $x$

1:$i=0$

2:Compute $r_{0}=b-Ax_{0}$ $\triangleright$ $u$

3:while not converged do

4: Compute $\begin{bmatrix}f_{i}\\\ g_{i}\end{bmatrix}=\begin{bmatrix}b\\\
0\end{bmatrix}-\begin{bmatrix}I&A\\\
A^{T}&0\end{bmatrix}\begin{bmatrix}r_{i}\\\ x_{i}\end{bmatrix}$
$\triangleright$ $u_{r}$

5: Solve $\min_{\delta x_{i}^{(1)}}\|f_{i}-A\delta x_{i}^{(1)}\|_{2}$
$\triangleright$ $u$

6: Solve $\min_{\delta x_{i}^{(2)}}\|r_{i}-A\delta x_{i}^{(2)}\|_{2}$
$\triangleright$ $u$

7: Solve $\min_{\delta r_{i}^{(1)}}\|g_{i}-A^{T}\delta r_{i}^{(1)}\|_{2}$
$\triangleright$ $u$

8: Compute $\delta r_{i}^{(2)}=P_{N}f_{i}$ $\triangleright$ $u$

9: Update $\begin{bmatrix}r_{i+1}\\\
x_{i+1}\end{bmatrix}=\begin{bmatrix}r_{i}\\\
x_{i}\end{bmatrix}+\begin{bmatrix}\delta r_{i}^{(1)}+\delta r_{i}^{(2)}\\\
\delta x_{i}^{(1)}+\delta x_{i}^{(2)}\end{bmatrix}$ $\triangleright$ $u$

10: $i=i+1$

11:end while

##  7 Numerical experiments

In order to compare the various methods numerically and to confirm the
behavior described by the analysis, we consider a simple class of least-
squares problems using MATLAB R2021a***Our code is available at
https://github.com/dauzickaite/LSIR. A $10^{3}\times 10$ coefficient matrix
$A$ is generated as

$A=\text{gallery}(\text{'randsvd'},[1e3,1e1],k,3)$

and we change the $\kappa(A)$ values (the $k$ parameter above). We generate
the right-hand side $b$ as

$b=mp(A)mp(y)+mp(e),$

where $mp(.)$ simulates precision accurate to 64 decimal digits via the
Advanpix toolbox [1], $y$ is a normalized random vector with uniformly
distributed entries, and $e$ is such that $A^{T}e=0$. We then cast $b$ to
double or single. We vary the size of $e$ and thus obtain results for
different residual norms. The true solution $x^{*}$ and true residual $r^{*}$
are computed as

$\displaystyle x^{*}=$ $\displaystyle\,mp(A)\backslash mp(b);$ $\displaystyle
r^{*}=$ $\displaystyle\,mp(b)-mp(A)*x^{*}.$

Note that because we are aiming for the unit roundoff level of accuracy, the
results are very sensitive to the setup. We explore two sets of precisions
$u_{r}$ for computing the residuals and $u$ for other computations: $u_{r}$
set to quadruple (simulated via Advanpix) and $u$ set to double, and $u_{r}$
set to double and $u$ set to single.

We declare that LSIR converged if

$\|x_{i}-x^{*}\|/\|x^{*}\|\leq\tau.$

We set $\tau=8u$, where $u$ is the unit roundoff for double or single
precision. If this condition is satisfied with $x_{0}$, then LSIR is not
performed. LSIR is run for a maximum of 30 iterations.

###  7.1 Solving via QR

We compute the full QR decomposition in the working precision $u$ as
$[\widetilde{Q},\widetilde{R}]=qr(A)$ and set $R=\widetilde{R}(1:n,1:n)$ and
$Q_{1}=\widetilde{Q}(:,1:n)$. The initial estimate is
$x_{0}=R\backslash(Q_{1}^{\prime}*b)$, where $Q_{1}^{\prime}$ denotes
$Q_{1}^{T}$. We show when the condition for recognizing the true solution in
the LS approach is satisfied in Figure 1. We use condition (37) to estimate
for which instances the augmented system LSIR approach is expected to
converge; see Figure 2. Note that we expect the semi-normal equations approach
to be insensitive to $\|r\|_{2}$ and converge for problems with
$\kappa(A)<u^{-1/2}$.

Fig. 1: LS system approach. Criteria for the least-squares system approach to
recognize a true solution if it is found (3) for a set of LS problems with
varying $\kappa(A)$ and $\|r^{*}\|$. The white cells indicate that (3) is
satisfied.

(a) $u_{r}$ quad, $u$ double

(b) $u_{r}$ double, $u$ single

Fig. 2: Augmented system approach. LSIR convergence criteria (37) for a set of
LS problems with varying $\kappa(A)$ and $\|r^{*}\|$ for two sets of
precisions $u_{r}$ and $u$. The white cells indicate that (37) is satisfied.

In Figures 3, 4, and 5, we show the iteration count, and relative errors
$\|x_{i}-x^{*}\|/\|x^{*}\|$ and $\|r_{i}-r^{*}\|/\|r^{*}\|$, respectively, for
the LS, semi-normal equations, and augmented system approaches. The augmented
system approach converges in all the cases where the criteria (37) is
satisfied (as shown in Figure 2) and even for more ill-conditioned systems
with larger residuals. As suggested by the theoretical analysis, the semi-
normal equations approach converges for all problems when $\kappa(A)<u^{-1/2}$
(and even for some worse conditioned ones) and is not sensitive to the size of
the residual.

Comparing all three approaches, the augmented system approach is able to
refine the solution $x$ for both more ill-conditioned problems and problems
with larger true residual. The semi-normal equations approach outperforms the
LS approach but not the augmented system approach. Note that in most cases
very few LSIR iterations are needed; in general a larger number of iterations
is needed when $\kappa(A)$ grows so large that the theoretical results cannot
guarantee convergence.

Our stopping criteria for LSIR is based on the error in $x$ only and we thus
may stop the algorithm too early to refine $r$ to the attainable level.
However, in the cases where the augmented system approach is unable to refine
the solution $x$, the error in the relative residual reaches the
$\mathcal{O}(u)$ level. The LS system approach mostly fails to refine the
residual to such accuracy even when $\|x_{i}-x^{*}\|/\|x^{*}\|$ reaches
$\mathcal{O}(u)$. The semi-normal equations approach refines the residual for
problems with small $\|r^{*}\|$. We can clearly observe that $r$ is harder to
refine with all the methods when $\|r^{*}\|$ decreases.

(a) LS approach, (quad, double)

(b) LS approach, (double, single)

(c) semi-normal, (quad, double)

(d) semi-normal, (double, single)

(e) augmented approach, (quad, double)

(f) augmented approach, (double, single)

Fig. 3: Count of iterative refinement iterations for the LS, semi-normal
equations, and augmented system approach when $u_{r}$ is set to quadruple and
$u$ is set to double (left panels), and $u_{r}$ is set to double and $u$ is
set to single (right). Black denotes when LSIR does not converge in 30
iterations.

(a) LS approach, (quad, double)

(b) LS approach, (double, single)

(c) semi-normal, (quad, double)

(d) semi-normal, (double, single)

(e) augmented approach, (quad, double)

(f) augmented approach, (double, single)

Fig. 4: As in Figure 3, but for the relative error in $x$.

(a) LS approach, (quad, double)

(b) LS approach, (double, single)

(c) semi-normal, (quad, double)

(d) semi-normal, (double, single)

(e) augmented approach, (quad, double)

(f) augmented approach, (double, single)

Fig. 5: As in Figure 3, but for the relative error in $r$.

###  7.2 Solving via iterative methods

We perform some numerical experiments solving the correction equations via
iterative methods, namely, LSQR for the least-squares system approach and
GMRES for the augmented system approach. We also test the combined approach in
Algorithm 6.1 using LSQR for the least-squares problems and setting
$P_{N}=I-Q_{1}Q_{1}^{T}$, where $Q_{1}$ is as in the previous section. Note
that the success of these iterative methods depends on factors such as
preconditioning, convergence criteria, the chosen scaling factor for the
augmented system, precisions for computing matrix-vector products and applying
the preconditioner, etc. We do not aim to provide an extensive study of these
effects but rather to illustrate some of the behavior. We thus perform
experiments for two cases: without preconditioning, and when a high quality
preconditioner is employed.

We use the built-in MATLAB implementation of LSQR and our own implementation
of GMRES. The convergence tolerances are set to $10^{-14}$ and $10^{-7}$ for
LSQR and $10^{-12}$ and $10^{-6}$ for GMRES when $u$ is set to double and
single, respectively. Additional experiments show that larger tolerances lead
to worse performance (results not shown). The maximum number of LSQR
iterations is set to $n=10$, and $50$ for GMRES. The initial approximate
solution $x_{0}$ is obtained via an unpreconditioned LSQR solve in all cases
and is then used to compute $r_{0}=b-Ax_{0}$. As was observed by Golub when
using a direct solver [15], additional numerical experiments show that the LS
system approach shows large sensitivity to the quality of $x_{0}$ (results not
shown here).

To generate the preconditioners, we compute the economic QR factorization of
$A$ as in the previous section. The $R$ factor is used as a right
preconditioner in LSQR for all the least-squares problems, except step 7 in
Algorithm 6.1 where it is used on the left. The augmented system is
preconditioned on the left as in [11], that is, we solve

$\begin{pmatrix}I-Q_{1}Q_{1}^{T}&Q_{1}R^{-T}\\\
R^{-1}Q_{1}^{T}&-R^{-1}R^{-T}\end{pmatrix}\begin{pmatrix}I&A\\\
A^{T}&0\end{pmatrix}\begin{pmatrix}\delta r_{i}\\\ \delta
x_{i}\end{pmatrix}=\begin{pmatrix}I-Q_{1}Q_{1}^{T}&Q_{1}R^{-T}\\\
R^{-1}Q_{1}^{T}&-R^{-1}R^{-T}\end{pmatrix}\begin{pmatrix}f_{i}\\\
g_{i}\end{pmatrix}.$

Results for unpreconditioned solves in Figures 6, 7, 8, and 9 show that all
the LSIR approaches except the augmented approach perform poorly. The LS and
combined approaches converge only for well conditioned problems. The augmented
approach converges in a few IR iterations for all problem where
$\kappa(A)^{2}$ is safely less than $u^{-1}$, although at a high cost of many
GMRES iterations; note that in our example $\sigma_{min}(A)=\kappa(A)^{-1}\ll
1$ and thus we can expect $\kappa(\widetilde{A})\approx\kappa(A)^{2}$ and the
results agree with the IR convergence theory for linear systems of equations.

Preconditioning for the iterative methods improves the situation
significantly; see Figures 10, 11, 12, and 13. The LS and augmented approaches
converge for similar sets of problems as when solving via QR decomposition.
The preconditioned combined approach performs similar to the preconditioned LS
approach, except the combined approach gives a substantially better
improvement for the residual at the cost of more LSQR iterations.

(a) LS approach, (quad, double)

(b) LS approach, (double, single)

(c) augmented approach, (quad, double)

(d) augmented approach, (double, single)

(e) combined approach, (quad, double)

(f) combined approach, (double, single)

Fig. 6: Count of iterative refinement iterations using unpreconditioned LSQR
and GMRES for the LS, augmented system, and combined approaches when $u_{r}$
is set to quadruple and $u$ is set to double (left panels), and $u_{r}$ is set
to double and $u$ is set to single (right). Black denotes when LSIR does not
converge in 30 iterations.

(a) LS approach, (quad, double)

(b) LS approach, (double, single)

(c) augmented approach, (quad, double)

(d) augmented approach, (double, single)

(e) combined approach, (quad, double)

(f) combined approach, (double, single)

Fig. 7: As in Figure 6, but for the cumulative number of inner iterative
refinement iterations, that is, LSQR iterations for the LS and combined
approaches, and GMRES iterations for the augmented approach. Black denotes
when LSIR does not converge in 30 iterations.

(a) LS approach, (quad, double)

(b) LS approach, (double, single)

(c) augmented approach, (quad, double)

(d) augmented approach, (double, single)

(e) combined approach, (quad, double)

(f) combined approach, (double, single)

Fig. 8: As in Figure 6, but for the relative error in $x$.

(a) LS approach, (quad, double)

(b) LS approach, (double, single)

(c) augmented approach, (quad, double)

(d) augmented approach, (double, single)

(e) combined approach, (quad, double)

(f) combined approach, (double, single)

Fig. 9: As in Figure 6, but for the relative error in $r$.

(a) LS approach, (quad, double)

(b) LS approach, (double, single)

(c) augmented approach, (quad, double)

(d) augmented approach, (double, single)

(e) combined approach, (quad, double)

(f) combined approach, (double, single)

Fig. 10: As in Figure 6, but when using preconditioned LSQR and GMRES.

(a) LS approach, (quad, double)

(b) LS approach, (double, single)

(c) augmented approach, (quad, double)

(d) augmented approach, (double, single)

(e) combined approach, (quad, double)

(f) combined approach, (double, single)

Fig. 11: As in Figure 10, but for the cumulative iterations of LSQR and GMRES.

(a) LS approach, (quad, double)

(b) LS approach, (double, single)

(c) augmented approach, (quad, double)

(d) augmented approach, (double, single)

(e) combined approach, (quad, double)

(f) combined approach, (double, single)

Fig. 12: As in Figure 10, but for the relative error in $x$.

(a) LS approach, (quad, double)

(b) LS approach, (double, single)

(c) augmented approach, (quad, double)

(d) augmented approach, (double, single)

(e) combined approach, (quad, double)

(f) combined approach, (double, single)

Fig. 13: As in Figure 10, but for the relative error in $r$.

##  8 Summary

We discussed the least-squares system, semi-normal equations, and augmented
system approaches to two-precision iterative refinement for least-squares
problems. The focus was on refining the solution $x$, although the augmented
system approach also explicitly refines the least-squares residual $r$. The
methods require a similar amount of computation in the high precision $u_{r}$;
see Table 2.

In addition to reformulating existing results to enable a comparison of the
different approaches, we provided a new analysis for the forward error in the
two-precision semi-normal equations approach. We further showed how to combine
the LS approach, which can use iterative least-squares solvers, and the
augmented system approach, which also refines the residual.

LSIR approach | $u_{r}$ | $u$
---|---|---
$Az$ | $A^{T}z$ | flops | flops
LS system | yes | no | $2mn$ | $4mn+n^{2}$
Semi-normal equations | yes | yes | $4mn-n$ | $2n^{2}+n$
Augmented system | yes | yes | $4mn+m-n$ | $8mn+2n^{2}$
Table 2: Computational cost per LSIR iteration for the three LSIR approaches
in two precisions $u_{r}$ and $u$ with $u_{r}\ll u$ when $A$ is dense and
Householder QR is used to solve the update equations. We specify the number of
flops and whether a matrix-vector product with $A$ and $A^{T}$ has to be
computed in the higher precision.

The main observations are summarized in the following.

  * •

We do not have theoretical convergence guarantees for the LS approach and know
that for some right-hand side vectors the approach may not be able to refine
the solution to the level of the working precision, namely, when the least-
squares residual is large. Numerical experiments show that the method does not
converge when $\kappa(A)>u^{-1/2}$.

  * •

The semi-normal equations approach converges for problems with $\kappa(A)\leq
u^{-1/2}$ and is not sensitive to the size of the least-squares residual.

  * •

Conditions for LSIR convergence in the solution and the residual when using
the augmented system approach are available. Convergence for problems with
high condition number (even when $\kappa(A)>u^{-1/2}$) and large residual can
be obtained with a direct solver; however, a suitable preconditioner is
crucial if iterative solvers are to be used. In either case, the situation can
be further improved by the use of an optimal scaling.

Thus the practitioner may choose to use the LS system approach if the problem
is well-conditioned and the least-squares residual is small. If it is known
that the problem is well-conditioned but the residual is large or unknown,
then the semi-normal equations approach may be used. If the user has no
knowledge of the conditioning and the size of the residual or it is known that
the problem is ill-conditioned, then the augmented system approach is the most
promising. The combined approach may be interesting in the cases where both
$x$ and $r$ are to be refined and parallel iterative least-squares solvers are
preferred.

## References

  * [1] Advanpix multiprecision computing toolbox for MATLAB.  http://www.advanpix.com.
  * [2] P. Amestoy, A. Buttari, N. J. Higham, J.-Y. L’Excellent, T. Mary, and B. Vieublé, Five-precision GMRES-based iterative refinement, SIAM Journal on Matrix Analysis and Applications, 45 (2024), pp. 529–552.
  * [3] M. Arioli, I. S. Duff, and P. P. de Rijk, On the augmented system approach to sparse least-squares problems, Numerische Mathematik, 55 (1989), pp. 667–684.
  * [4] Å. Björck, Iterative refinement of linear least squares solutions I, BIT Numerical Mathematics, 7 (1967), pp. 257–278.
  * [5] , Comment on the iterative refinement of least-squares solutions, Journal of the American Statistical Association, 73 (1978), pp. 161–166.
  * [6] , Stability analysis of the method of seminormal equations for linear least squares problems, Linear Algebra and its Applications, 88 (1987), pp. 31–48.
  * [7] , Numerical methods for least squares problems, SIAM, 1996.
  * [8] Å. Björck and C. C. Paige, Solution of augmented linear systems using orthogonal factorizations, BIT Numerical Mathematics, 34 (1994), pp. 1–24.
  * [9] P. Businger and G. H. Golub, Linear least squares solutions by Householder transformations, Numerische Mathematik, 7 (1965), pp. 269–276.
  * [10] E. Carson and N. J. Higham, Accelerating the solution of linear systems by iterative refinement in three precisions, SIAM Journal on Scientific Computing, 40 (2018), pp. A817–A847.
  * [11] E. Carson, N. J. Higham, and S. Pranesh, Three-precision GMRES-based iterative refinement for least squares problems, SIAM Journal on Scientific Computing, 42 (2020), pp. A4063–A4083.
  * [12] Y. Chen, X. Ji, and J. Nagy, Half-precision Kronecker product SVD preconditioner for structured inverse problems, arXiv preprint arXiv:2311.15393, (2023).
  * [13] J. Demmel, Y. Hida, E. J. Riedy, and X. S. Li, Extra-precise iterative refinement for overdetermined least squares problems, ACM Transactions on Mathematical Software (TOMS), 35 (2009), pp. 1–32.
  * [14] E. N. Epperly, Fast and forward stable randomized algorithms for linear least-squares problems, arXiv preprint arXiv:2311.04362, (2023).
  * [15] G. Golub, Numerical methods for solving linear least squares problems, Numerische Mathematik, 7 (1965), pp. 206–216.
  * [16] G. H. Golub and J. H. Wilkinson, Note on the iterative refinement of least squares solution, Numerische Mathematik, 9 (1966), pp. 139–148.
  * [17] N. J. Higham, Accuracy and Stability of Numerical Algorithms, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, second ed., 2002.
  * [18] N. J. Higham and S. Pranesh, Exploiting lower precision arithmetic in solving symmetric positive definite linear systems and least squares problems, SIAM Journal on Scientific Computing, 43 (2021), pp. A258–A277.
  * [19] H. Li, Double precision is not necessary for LSQR for solving discrete linear ill-posed problems, Journal of Scientific Computing, 98 (2024), p. 55.
  * [20] C. C. Paige and M. A. Saunders, LSQR: An algorithm for sparse linear equations and sparse least squares, ACM Transactions on Mathematical Software (TOMS), 8 (1982), pp. 43–71.
  * [21] Y. Saad, Iterative methods for sparse linear systems, SIAM, 2003.
  * [22] J. Scott and M. Tůma, A computational study of using black-box QR solvers for large-scale sparse-dense linear least squares problems, ACM Transactions on Mathematical Software (TOMS), 48 (2022), pp. 1–24.