# Are sketch-and-precondition least squares solvers numerically stable?

Maike Meier Mathematical Institute, University of Oxford.
<EMAIL_ADDRESS>Yuji Nakatsukasa Mathematical Institute, University of
Oxford<EMAIL_ADDRESS>Alex Townsend Department of Mathematics,
Cornell University<EMAIL_ADDRESS>Marcus Webb Department of
Mathematics, University of Manchester<EMAIL_ADDRESS>
###### Abstract

Sketch-and-precondition techniques are popular for solving large least squares
(LS) problems of the form $Ax=b$ with $A\in\mathbb{R}^{m\times n}$ and $m\gg
n$. This is where $A$ is “sketched” to a smaller matrix $SA$ with
$S\in\mathbb{R}^{\lceil cn\rceil\times m}$ for some constant $c>1$ before an
iterative LS solver computes the solution to $Ax=b$ with a right
preconditioner $P$, where $P$ is constructed from $SA$. Popular sketch-and-
precondition LS solvers are Blendenpik and LSRN. We show that the sketch-and-
precondition technique is not numerically stable for ill-conditioned LS
problems. Instead, we propose using an unpreconditioned iterative LS solver on
$(AP)y=b$ with $x=Py$ when accuracy is a concern. Provided the condition
number of $A$ is smaller than the reciprocal of the unit round-off, we show
that this modification ensures that the computed solution has a comparable
backward error to the iterative LS solver applied to a well-conditioned
matrix. Using smoothed analysis, we model floating-point rounding errors to
provide a convincing argument that our modification is expected to compute a
backward stable solution even for arbitrarily ill-conditioned LS problems.

###### keywords:

Least squares, numerical stability, sketching, preconditioner

65F10, 65F20

##  1 Introduction

Randomized numerical linear algebra is a growing subfield of matrix
computations that has produced major advances in low-rank approximation [14],
iterative methods [21], and projections [4]. Sketch-and-precondition
techniques are a class of randomized algorithms for solving overdetermined
least squares (LS) problems of the form

(1) $\min_{x\in\mathbb{R}^{n}}\|Ax-b\|_{2},\qquad A\in\mathbb{R}^{m\times
n},\quad b\in\mathbb{R}^{m\times 1},$

where $m>n$. One first sketches $A$ to a smaller matrix $SA$ with
$S\in\mathbb{R}^{\lceil cn\rceil\times m}$ for some constant $c>1$, then a
right preconditioner, $P$, is constructed from $SA$. Finally, one solves
$Ax=b$ using an iterative LS solver with the right preconditioner $P$. There
are many details of sketch-and-precondition based on how to sketch and
construct $P$ as well as which iterative LS solver to employ. One of the most
prominent sketch-and-precondition techniques is known as Blendenpik [1] (see
Algorithm 1). In exact arithmetic, Blendenpik has a complexity of
$\mathcal{O}(mn\log m)$ operations, which is better than the
$\mathcal{O}(mn^{2})$ QR-based direct solver. Consequently, for large LS
problems, Blendenpik can be substantially faster than the LS solver
implemented in LAPACK [1]. However, are sketch-and-precondition
techniques—such as Blendenpik—numerically stable?

Surprisingly, we find that sketch-and-precondition techniques such as
Blendenpik and LSRN [16] are numerically unstable in their standard form (see
Fig. 1). However, a modification (see Algorithm 2) can ensure that the
computed residual is close to optimal and the backward error is approximately
machine precision. The modification is simple. Instead of using $P$ as a
preconditioner, we explicitly compute $AP$ and then employ an unpreconditioned
iterative LS solver on $(AP)y=b$ with $x=Py$. We call this a sketch-and-apply
technique. Of course, in exact arithmetic sketch-and-precondition and sketch-
and-apply compute the same solution; however, for ill-conditioned LS problems
in floating-point arithmetic, we find a significant difference. To our
knowledge, this is the first algorithm for LS problems based on randomized
sketching that is demonstrated to be backward stable (with a mild assumption
to be made precise in Theorem 3.9).

Unfortunately, while sketch-and-precondition techniques cost
$\mathcal{O}(mn\log m)$ operations, sketch-and-apply costs
$\mathcal{O}(mn^{2})$ operations as $AP$ must be computed. Therefore, sketch-
and-apply techniques have the same computational complexity as the classical
QR-based LS solver. However, we show that they are competitive for large, ill-
conditioned LS problems (see Fig. 4). We also note that LSRN, another popular
sketch-and-precondition algorithm, also uses $\mathcal{O}(mn^{2})$ operations
but is still competitive as the expensive computation is in matrix
multiplication, which is highly parallelizable. The same applies to sketch-
and-apply.

A popular iterative algorithm for solving Eq. 1 is LSQR [18], which is
competitive when $A$ is large, sparse, and well-conditioned. The number of
iterations required to reach a desired accuracy typically depends on the
condition number of $A$ [2, Ch. 7.4]. As a result, a preconditioner is
necessary when $A$ is ill-conditioned to ensure a reasonable speed of
convergence. With careful sketching, the condition number of $AP$ with $P$
from Algorithm 1 is small with high probability in exact arithmetic [19, Lem.
1]. In finite precision arithmetic, the condition number of the computed $AP$
remains modest in size, provided that $\kappa_{2}(A)u\ll 1$, where
$\kappa_{2}(A)$ is the condition number of $A$ and $u$ is the unit round-off
(see Theorem 3.7).

Algorithm 1 A sketch-and-precondition LS solver for Eq. 1 known as Blendenpik.
Here, HHQR refers to the Householder QR algorithm.

1: Draw a random sketching matrix $S\in\mathbb{R}^{s\times m}$, where $m\gg
s>n$

2: Compute $B=SA$

3: Compute the triangular factor $R$ of a QR factorization of $B$ using HHQR

4: Solve $Ax=b$ with LSQR and right preconditioner $P=R^{-1}$

However, when $A$ is an extremely ill-conditioned matrix, so is $P=R^{-1}$,
and sketch-and-precondition solves linear systems with an ill-conditioned $R$,
leading to a poor final LS residual in finite precision. Indeed, it is
recommended in the Blendenpik paper that LSQR is avoided when
$\kappa_{2}(R)>1/(5u)$ [1]. What we find is that for moderately ill-
conditioned systems (such as $\kappa_{2}(A)\approx\sqrt{u^{-1}}$), the
residual of the computed solution is affected commensurately (see Section
1.1).

However, in sketch-and-apply, we explicitly compute $AP$ before using an
unpreconditioned LSQR. We show that this leads to a final computed LS residual
that is significantly more accurate (see Algorithm 2), especially when
$\kappa_{2}(A)$ is extremely large. In particular, we prove that if
$\kappa_{2}(A)\ll u^{-1}$, our sketch-and-apply technique computes a backward
stable solution provided that LSQR on a well-conditioned matrix computes a
backward stable solution (see Theorem 3.9).

Algorithm 2 A sketch-and-apply LS solver for Eq. 1. Here, HHQR refers to the
Householder QR algorithm.

1: Draw a random sketching matrix $S\in\mathbb{R}^{s\times m}$, where $m\gg
s>n$

2: Compute $B=SA$

3: Compute the triangular factor $R$ of a QR factorization of $B$ using HHQR

4: Compute $Y=AR^{-1}$ with forward substitution

5: Solve $Yz=b$ with LSQR and no preconditioner

6: Compute $x=R^{-1}z$ with back substitution

When $\kappa_{2}(A)\gtrsim u^{-1}$, our analysis does not guarantee that
Algorithm 2 leads to an accurate solution. However, in practice, we often get
accurate final LS residuals. To explain why this is the case, we model
floating-point rounding errors using smoothed analysis. That is, we consider
“smoothing” $A$ to $A+\sigma G/\sqrt{m}$, where the entries of $G$ are
independent and identically distributed (i.i.d.) standard Gaussian random
variables and $\sigma$ is a scaling factor that we select as
$\sigma=10\|A\|_{2}u$. The idea is that $A+\sigma G/\sqrt{m}$ is significantly
better conditioned than $A$ itself, assuming that $A$ is extremely ill-
conditioned. In fact, for sufficiently small $\sigma$, one can show that
$\kappa_{2}(A+\sigma G/\sqrt{m})\lesssim 1/\sigma$ with high probability [5].
The additive perturbation of $\sigma G/\sqrt{m}$ to $A$ ensures that sketch-
and-apply techniques can also compute solutions with good backward error, even
when $A$ is extremely ill-conditioned (see Algorithm 3). For extremely ill-
conditioned LS problems, one could explicitly add an additive random
perturbation or hope that floating point rounding errors deliver the same
effect, as it often does in practice.

Algorithm 3 A smoothed sketch-and-apply LS solver for Eq. 1 when
$\kappa_{2}(A)\gtrsim u^{-1}$.

1: Draw a random standard Gaussian matrix $G\in\mathbb{R}^{m\times n}$ with
i.i.d. entries

2: Compute $\tilde{A}=A+\sigma G/\sqrt{m}$ for $\sigma=10\|A\|_{2}u$

3: Perform Algorithm 2 on $\tilde{A}$

Throughout the paper, we assume that $A\in\mathbb{R}^{m\times n}$ and
$SA\in\mathbb{R}^{s\times n}$, $m>s>n$, are of full rank, so that Eq. 1 has a
unique solution.

###  1.1 The numerical instabilities of sketch-and-precondition

We begin by demonstrating the numerical instabilities that occur with sketch-
and-precondition techniques.111Experiments are performed in MATLAB 2022b using
64-bit arithmetic on a single core of a MacBook Pro equipped with a 2.3 GHz
Dual-Core Intel Core i5 processor and 16 GB 2133 MHz LPDDR3 of system memory.
We construct LS problems at random by setting $A=U\Sigma V^{T}$, where
$U\in\mathbb{R}^{m\times m}$ and $V\in\mathbb{R}^{n\times n}$ are random
orthogonal matrices from a Haar distribution and
$\Sigma=\mbox{diag}(\sigma_{1},\ldots,\sigma_{n})$, where
$\sigma_{1},\ldots,\sigma_{n}$ are logarithmically spaced between $1$ and
$\xi$. We set $m=2000$, $n=100$, and take three different values of
$\xi=0.1,10^{-10},\text{ and, }10^{-15}$. The right-hand side of the LS
problem is generated as $b=Ax^{*}+e$, where $x^{*}$ is a random vector with
i.i.d. standard Gaussian entries and $e$ represents a noise vector. The noise
is chosen to be orthogonal to the column space of $A$ so that $x^{*}$ is the
exact solution to $Ax=b$. For this reason, we do not expect the computed
residual to go below $\|e\|_{2}$. For each value of $\xi$, we consider two
values of $\|e\|_{2}$: (i) $\|e\|_{2}=10^{-3}$ (see Fig. 1 (left column)) and
(ii) $\|e\|_{2}=10^{-14}$ (see Fig. 1 (right column)).

To solve the LS problems, we use (1) sketch-and-precondition Blendenpik (see
Algorithm 1) and (2) sketch-and-apply Blendenpik (see Algorithm 2). The only
difference between the two algorithms is how they use the preconditioner. We
use a sketching matrix $S\in\mathbb{R}^{8n\times m}$, known as a subsampled
randomized cosine transform(see Section 2.2.2).

Figure 1: Sketch-and-precondition Blendenpik (circles) and sketch-and-apply
Blendenpik (crosses) applied to our noisy and randomly generated LS problems
where $\kappa_{2}(A)=10$ (blue), $\kappa_{2}(A)=10^{10}$ (red), and
$\kappa_{2}(A)=10^{15}$ (black). We take two different noise levels of
$\|e\|_{2}=10^{-3}$ (left column) and $\|e\|_{2}=10^{-14}$ (right column).

Fig. 1 compares the performance of both algorithms in terms of the relative
residual $\|Ax-b\|_{2}/\|b\|_{2}$ (top row), the relative forward error
$\|x-x^{*}\|_{2}/\|x^{*}\|_{2}$ (middle row), and the normwise backward error
$\eta_{F}(x)$ (bottom row, see Eq. 2). We consider these three quantities as
each one is an important performance measure. A stable algorithm should
compute a final solution with a relative residual close to the noise level
given by $\|e\|_{2}$, a forward error approximately equal to
$\kappa_{2}(A)u(1+\kappa_{2}(A)\|e\|_{2})$ [11, 25], and a backward error of
order $u$, where the backward error is given by [23] (also, see [11, Thm.
20.5])

(2) $\eta_{F}(x):=\min\\{\|[\Delta A,\Delta b]\|_{F}\,:\,\|(A+\Delta
A)x-(b+\Delta b)\|_{2}=\min\\}.$

Based on the final residual, forward error, and backward error, we find that
sketch-and-precondition Blendenpik underperforms on LS problems with matrices
$A$ that are moderately or highly ill-conditioned; especially when the noise
level is low.

To examine things more closely, consider the final relative residuals (top row
of Fig. 1). When the noise level is small, sketch-and-apply always converges
rapidly to the solution, regardless of $\kappa_{2}(A)$. However, sketch-and-
precondition does not converge to the solution, and the final accuracy seems
proportional to $\kappa_{2}(A)u$. This is a floating-point rounding effect
because both algorithms employ the same preconditioner in different ways. We
see a similar story for the relative forward error (middle row of Fig. 1) with
ill-conditioned LS problems, causing sketch-and-precondition Blendenpik to
stagnate before the desired accuracy.

Lastly, and most importantly for numerical analysts, sketch-and-precondition
Blendenpik does not compute a solution with a backward error that is close to
machine precision. In contrast, sketch-and-apply does (bottom row of Fig. 1).
We conclude that sketch-and-precondition Blendenpik is not numerically stable.

Figure 2: A comparison of the convergence of the relative residual for sketch-
and-precondition Blendenpik, sketch-and-apply Blendenpik, unpreconditioned
LSQR, and preconditioned LSQR with $R^{-1}$ from the QR factorization of $A$.

The numerical instabilities in sketch-and-precondition Blendenpik could be due
to multiple sources. To demonstrate that it is the way that the preconditioner
is employed, and not how it is formed, we try one more randomly generated LS
problem with $A\in\mathbb{R}^{2000\times 100}$, $\kappa_{2}(A)=10^{12}$, and a
noise level of $\|e\|_{2}=10^{-14}$. We solve the LS problem with the
unpreconditioned LSQR and preconditioned LSQR, where the preconditioner is
obtained by computing the QR factorization of $A$. The preconditioned LSQR
approach is not practical because it requires the QR factorization of $A$. We
still find in Fig. 2 that sketch-and-apply Blendenpik is the only approach
that computes a solution with a residual close to machine precision.

While we have not spotted the precise source of instability, we believe it
lies in the fact that each application of the preconditioner $P$—which
involves solving an ill-conditioned linear system $Rx=b$—incurs a relative
error proportional to $u\kappa_{2}(R)$. Such errors are present every time we
apply $P$ or $P^{T}$. Moreover, each application behaves somewhat differently,
in that in the $i$th iterate we have $(R+\Delta R_{i})\hat{x}_{i}=b_{i}+\delta
b_{i}$ where $\Delta R_{i},\delta b_{i}$ are small but different for each $i$.
In other words, one can view the preconditioner as having an $u\kappa_{2}(R)$
nonlinear effect. The fact that the convergence of iterative methods can get
impaired by nonlinear preconditioners has been observed in [24].

###  1.2 Paper structure

In Section 2, we introduce some background material on sketch-and-precondition
LS solvers, numerical stability analysis, and smoothed analysis. In Section 3,
we consider sketch-and-apply Blendenpik (see Algorithm 2) in finite precision.
In Section 4, we consider extremely ill-conditioned LS problems and look at
the numerical stability of smoothed sketch-and-apply Blendenpik (see Algorithm
3). We introduce a master algorithm and display numerical experiments in
Section 5. Finally, in Section 6, we argue that despite the
$\mathcal{O}(mn^{2})$ complexity of sketch-and-apply, the technique is
competitive.

While in this paper we focus on LS problems, we expect much of the stability
results to carry over to solving underdetermined linear systems using
sketching, as done in LSRN [16].

##  2 Background material

We now introduce some background material for sketch-and-precondition
techniques (see Section 2.1), random sketching matrices (see Section 2.2),
numerical stability analysis (see Section 2.3), and smoothed analysis (see
Section 2.4).

###  2.1 Sketch-and-precondition least squares solvers

The main idea behind the sketch-and-precondition technique is that a
reasonable preconditioner for the LS problem in Eq. 1 can be constructed from
a sketch of $A$. In Blendenpik, the matrix $A$ is sketched to a small tall-
skinny matrix $SA$ and the preconditioner is taken to be the inverse of the
upper triangular factor from a QR factorization of $SA$, i.e., $P=R^{-1}$. In
exact arithmetic, the condition number of $AP$ is equal to $SQ_{A}$, where
$Q_{A}$ is an orthonormal basis for the column space of $A$.

######  Lemma 2.1.

Let $A\in\mathbb{R}^{m\times n}$ have linearly independent columns,
$S\in\mathbb{R}^{s\times m}$ with $s\geq n$ have linearly independent rows,
and $B=SA$. If $A=Q_{A}R_{A}$ and $B=QR$ are economized QR factorizations of
$A$ and $B$, respectively, then

$\kappa_{2}(AP)=\kappa_{2}(SQ_{A}).$

where $P=R^{-1}$.

######  Proof 2.2.

The proof follows the same argument by Meng, Saunders, and Mahoney in [16,
Lemma 4.2] and Rokhlin and Tygert in [19, Thm. 1] Note $AP=Q_{A}R_{A}R^{-1}$
so that $\kappa_{2}(AP)=\kappa_{2}(R_{A}R^{-1})$. Now, since
$B=SA=QR=S(Q_{A}R_{A})$, we have that $(SQ_{A})R_{A}R^{-1}=Q$ is orthonormal.
Hence, $\kappa_{2}(R_{A}R^{-1})=\kappa_{2}(SQ_{A})$.

Lemma 2.1 shows that the main idea behind sketch-and-precondition techniques
is excellent in exact arithmetic. In particular, the value of $\kappa_{2}(AP)$
is independent of $\kappa_{2}(A)$. In particular, one expects rapid
convergence of iterative LS solvers using the preconditioner $P$. While this
is correct in exact arithmetic, we have seen that rounding errors cause
significant problems for ill-conditioned LS problems (see Section 1.1).
Regardless, in exact arithmetic, we are left with the task of designing a
sketching matrix so that $\kappa_{2}(SQ_{A})$ is close to $1$ with high
probability.

###  2.2 Sketching matrices

Given Lemma 2.1, it is paramount to understand how to construct a sketching
matrix $S$ so that $SQ_{A}$ is well-conditioned. Matrices that achieve this
are called sketching matrices. In general, we desire $m\gg s>n$ so that $SA$
has full column rank and computing a preconditioner from $SA$ is
computationally efficient.

There are various ways to construct sketching matrices; almost all are
randomly generated. We consider two important types here: (1) Gaussian
matrices [16, 22, 26], and (2) subsampled randomized trigonometric transforms
(SRTTs) [1, 14, 19, 20].

####  2.2.1 Gaussian sketching matrices

A Gaussian sketching matrix of size $s\times m$ is a matrix with i.i.d.
Gaussian entries of mean $0$ and variance $1/s$. If one selects $s\geq\lceil
cn\rceil$ with $c>1$, then one can show that $\kappa_{2}(SQ_{A})$ is a small
constant depending on $c$ with high probability [14]. The drawback of Gaussian
matrices is the cost of computing $SA$. In particular, $SA$ costs
$\mathcal{O}(mn^{2})$ operations as $s=\mathcal{O}(n)$. However, Gaussian
sketching matrices are often used in theoretical analysis due to the
availability of excellent probability theory and concentration of measure
arguments [10, 14].

####  2.2.2 Subsampled randomized trigonometric transforms

An SRTT has the form $S=\mathcal{S}FD$, where $D$ is a square diagonal matrix
with $\pm 1$ entries at random, $F$ is an orthogonal trigonometric transform
(e.g., Fourier, cosine, or Hadamard), and $\mathcal{S}$ is an $s\times m$
scaled sampling matrix with one non-zero entry per row. For an SRTT, if
$s\geq\lceil cn\log n\rceil$ with some constant $c>1$, then $SQ_{A}$ is well-
conditioned with high probability [20]. In practice, the $\log n$ factor in
$s$ can often be dropped [14]. Moreover, the matrix-matrix product $SA$ can be
computed in $\mathcal{O}(mn\log s)$ operations, although it is often
computationally more efficient to compute $SA$ as $\mathcal{S}(F(DA))$
requiring $\mathcal{O}(mn\log m)$ operations. We use subsampled randomized
cosine transforms as the sketching matrices in our numerical experiments.

###  2.3 Numerical stability analysis

Basic arithmetic operations (e.g., $+$, $-$, $\cdot$, and $/$) on a computer
are performed with rounding errors due to the finite precision of numbers in
floating-point representation. The stability of a numerical algorithm refers
to the property to compute solutions with a small backward error [11] in the
presence of rounding errors. We are particularly interested in proving that
sketch-and-apply Blendenpik computes backward stable solutions. The backward
error for LS problems is defined in Eq. 2. A backward stable algorithm
computes solutions with a backward error of the same order as the precision.

We assume that computations are performed with a precision of $u\ll 1$ and
that numbers and arithmetic operations are exact up to this precision. We now
consider three example algorithms: (1) The standard algorithm for matrix-
matrix multiplication, (2) Householder QR, and (3) Triangular solve with
substitution. We need these results to understand sketch-and-apply LS solvers.

####  2.3.1 Matrix-matrix multiplication

Consider matrix-matrix multiplication of two matrices $A\in\mathbb{R}^{m\times
n}$ and $B\in\mathbb{R}^{n\times p}$. It is shown that the standard algorithm
for computing $AB$ satisfies [11, Chapt. 3]

(3) $|AB-{\rm
fl}(AB)|\leq\frac{nu}{1-nu}|A||B|=\gamma_{n}|A||B|,\qquad\gamma_{n}=\frac{nu}{1-nu},$

where ${\rm fl}(AB)$ means that $AB$ is computed with precision $u$ and $|A|$
denotes the entry-wise absolute value of $A$. The inequality in Eq. 3 is
understood as holding for each entry. It can be shown that Eq. 3 leads to the
bound

(4) $\|AB-{\rm
fl}(AB)\|_{2}\leq\gamma_{n}\min(\sqrt{m},\sqrt{n})\min(\sqrt{n},\sqrt{p})\|A\|_{2}\|B\|_{2},$

where $\|\cdot\|_{2}$ is the spectral norm. The inequality in Eq. 4 is a
forward error in the sense that it shows that ${\rm fl}(AB)$ is close to $AB$.

####  2.3.2 Householder QR factorization

Next, we consider the Householder QR of a matrix $A\in\mathbb{R}^{m\times n}$,
where $m\geq n$. We denote the triangular factor computed by the Householder
QR algorithm as $\hat{R}\in\mathbb{R}^{n\times n}$, and we are interested in
the accuracy of $\hat{R}$. By [11, Thm. 19.4], there exists a matrix
$Q\in\mathbb{R}^{m\times n}$ with orthonormal columns such that

$A+E=Q\hat{R},\qquad\|e_{j}\|_{2}\leq\tilde{\gamma}_{mn}\|a_{j}\|_{2},\quad
1\leq j\leq n,$

where $e_{j}$ and $a_{j}$ denote the $j$th columns of $E$ and $A$
respectively, and $\tilde{\gamma}_{mn}=cmnu/(1-cmnu)$ for a small integer
constant $c$. A direct consequence is that

(5) $\|E\|_{F}\leq\tilde{\gamma}_{mn}\|A\|_{F},$

where $\|\cdot\|_{F}$ is the matrix Frobenius norm. The inequality in Eq. 5
has the interpretation that the QR factorization computed by Householder QR is
exact for a slightly perturbed matrix $A$.

####  2.3.3 Triangular system solving with substitution

Finally, we consider the error arising from solving a triangular system with
substitution. Let $R\in\mathbb{R}^{n\times n}$ be a nonsingular upper-
triangular matrix. It is known that the computed solution, $\hat{x}$, of the
linear system $R^{T}x=b$ satisfies [11, Sec. 8.2]

(6)
$\|x-\hat{x}\|_{2}\leq\frac{\sqrt{n}\gamma_{n}\kappa_{2}(R)}{1-\sqrt{n}\gamma_{n}\kappa_{2}(R)}\|x\|_{2},$

where $\gamma_{n}$ is defined in Eq. 3. The inequality in Eq. 6 tells us that
we expect the relative forward error of the computed solution to be on the
order of $\kappa_{2}(R)$.

When we do any sketch-and-apply technique, we need to compute $AP$, where
$P=R^{-1}$. We compute $AP$ by solving $R^{T}x_{i}=a_{i}$ for $1\leq i\leq m$,
where $a_{i}^{T}$ denotes the $i$th row of $A$ and $x_{i}^{T}$ is the $i$th
row of $AP$. If we denote the computed solution to $R^{T}x_{i}=a_{i}$ by
$\hat{x}_{i}$, then from Eq. 6 we find that

$\displaystyle\sum_{i=1}^{m}\|x_{i}-\hat{x}_{i}\|_{2}^{2}\leq\left[\frac{\sqrt{n}\gamma_{n}\kappa_{2}(R)}{1-\sqrt{n}\gamma_{n}\kappa_{2}(R)}\right]^{2}\sum_{i=1}^{m}\|x_{i}\|_{2}^{2}.$

We conclude that for $P=R^{-1}$ we have

(7)
$\|AP-\widehat{AP}\|_{F}\leq\frac{\sqrt{n}\gamma_{n}\kappa_{2}(P)}{1-\sqrt{n}\gamma_{n}\kappa_{2}(P)}\|AP\|_{F},$

where $\widehat{AP}$ denotes the computed matrix-matrix product. Roughly
speaking, $\widehat{AP}$ is close to $AP$ when $\kappa_{2}(P)$ is modest. We
additionally bound the backward error resulting from solving $Rx=b$ with back
substitution. The computed solution $\hat{x}$ satisfies [11, Theorem 8.5]

(8)
$(\hat{R}+\Delta\hat{R})\hat{x}=\hat{z},\quad\|\Delta\hat{R}\|_{2}\leq\gamma_{n}\sqrt{n}\|\hat{R}\|_{2}.$

###  2.4 Smoothed analysis of condition numbers

Due to rounding errors, most matrices with $\kappa_{2}(A)>u^{-1}$ are
perturbed to a matrix with condition number less than $u^{-1}$ once
represented on a computer. To explain this, we can model the rounding errors
by an additive Gaussian perturbation to $A$. This type of technique fits into
the field of smoothed analysis and has been studied in [12]. More precisely,
we suppose that $A\in\mathbb{R}^{m\times n}$ and that we would like to bound
the condition number of $A+\sigma G$, where $\sigma$ is a small number and $G$
is a standard Gaussian matrix with i.i.d. entries. We have the following
statement, which is a corollary of [5, Thm. 1.1].

######  Corollary 2.3.

Let $A\in\mathbb{R}^{m\times n}$ with $m\geq 3n$ and $r\geq 10^{4}$. If $A$ is
perturbed to $\hat{A}=A+\sigma G/\sqrt{m}$ with $\sigma=8.25\|A\|_{2}/r$,
where $G$ is standard Gaussian matrix with i.i.d. entries, then

$\mathbb{P}\\!\left[\kappa_{2}(\hat{A})\geq r\,\right]<2^{n-m}.$

######  Proof 2.4.

The result immediately follows by setting $z=(r(1-\lambda))/e$ and
$\sigma=8.25/(r\sqrt{m})$ in Theorem 1.1 of [5]. This takes into account that
Theorem 1.1 of [5] considers the case when $m<n$. The extra $\|A\|_{2}$ factor
of $\sigma$ in the statement of the corollary ensures that our final result
continues to hold when $\|A\|_{2}\neq 1$.

Corollary 2.3 shows us that a small additive Gaussian perturbation to a matrix
$A$ ensures that the condition number is small with high probability. It
partially explains why matrices represented in floating-point arithmetic
rarely have a condition number $>u^{-1}$. In Section 4, we use Corollary 2.3
to explain why sketch-and-apply techniques continue to deliver accurate LS
solutions, even when we have $\kappa_{2}(A)>u^{-1}$.

##  3 On the numerical stability of sketch-and-apply Blendenpik

By the properties of sketching matrices and Lemma 2.1, we know that if one
computes $P$ and $AP$ in exact arithmetic, then $\kappa_{2}(AP)$ is bounded by
a constant independent of $\kappa_{2}(A)$ with high probability. In finite
precision, one typically needs to be more careful. In this section, we show
that when $\kappa_{2}(A)\ll u^{-1}$, the condition number of the computed
matrix $AP$ can also be bounded by a constant with high probability (see
Theorem 3.7). Unfortunately, such a bound does not extend to the case when
$\kappa_{2}(A)>u^{-1}$ (see Section 4).

There are three main steps to compute $AP$. In finite precision, the computed
matrices satisfy the following:

  1. 1.

$\hat{B}=SA+E_{1}$,

  2. 2.

$\tilde{Q}\hat{R}=\hat{B}+E_{2}$,

  3. 3.

$\hat{Y}=A\hat{R}^{-1}+\Delta Y$,

where hats denote that we are accounting for rounding errors. The matrix-
matrix multiplication error $E_{1}$ can be bounded using Eq. 4. We have a
bound on the error arising from the Householder QR algorithm $E_{2}$ in Eq. 5.
Note that $\tilde{Q}$ has orthonormal columns but is not equal to $Q$, where
$Q$ is the orthonormal factor of the exact QR factorization of $B=SA$.
Finally, the error term $\Delta Y$ in the step resulting from solving the
triangular matrix system can be bounded using Eq. 7. After all these errors,
we want to derive a bound on $\kappa_{2}(\hat{Y})$, as $\hat{Y}$ is the
computed preconditioned matrix in finite precision.

###  3.1 Bounding $\mathbf{\boldsymbol{\kappa}_{2}(\hat{Y})}$ in terms of
numerical errors

First, let $E=E_{1}+E_{2}$ be the sum of the additive errors from the matrix-
matrix multiplication and the QR factorization, such that
$\tilde{Q}\hat{R}=QR+E$. We show in Lemma 3.1 that the quantity
$\varepsilon_{1}:=\|E\|_{2}\|R^{-1}\|_{2}$ controls the numerical error in the
first two steps.

######  Lemma 3.1.

Let $A=Q_{A}R_{A}$ and $SA=QR$ be economized QR factorizations of $A$ and
$B=SA$, respectively. Let $\hat{B}$ be the computed $SA$ and $\hat{R}$ be the
computed upper triangular factor from the Householder QR algorithm of
$\hat{B}$, both computed in finite precision. Let $\tilde{Q}$ be a matrix with
orthonormal columns such that $E=\tilde{Q}\hat{R}-QR$ is the numerical error
associated with $\hat{R}$. Define $\varepsilon_{1}:=\|E\|_{2}\|R^{-1}\|_{2}$.
Assuming that $\varepsilon_{1}<1$, we find that $\hat{R}$ is invertible and
$A\hat{R}^{-1}$ is full rank. Furthermore,

$\kappa_{2}(\hat{R})\leq\frac{\kappa_{2}(R)+\varepsilon_{1}}{1-\varepsilon_{1}}\leq\frac{\kappa_{2}(SQ_{A})\kappa_{2}(A)+\varepsilon_{1}}{1-\varepsilon_{1}},$

and

$\kappa_{2}(A\hat{R}^{-1})\leq\kappa_{2}(SQ_{A})\frac{1+\varepsilon_{1}}{1-\varepsilon_{1}}.$

######  Proof 3.2.

We have

(9)
$\|\hat{R}\|_{2}=\|\tilde{Q}\hat{R}\|_{2}=\|QR+E\|_{2}\leq\|R\|_{2}+\|E\|_{2},$

and

(10)
$\sigma_{\min}(\hat{R})=\sigma_{\min}(QR+E)\geq\sigma_{\min}(R)-\|E\|_{2}.$

We see that $\sigma_{\min}(\hat{R})>0$ provided that $\varepsilon_{1}<1$. We
can combine Eq. 9 and Eq. 10 together to obtain an upper bound on
$\kappa_{2}(\hat{R})$. To bound $\kappa_{2}(A\hat{R}^{-1})$, we use the fact
that222The ‘$\dagger$’ superscript on a matrix denotes the matrix pseudo-
inverse. $A\hat{R}^{-1}\tilde{Q}^{T}=AR^{-1}(Q+ER^{-1})^{\dagger}$ to find
that for any $v\in\mathbb{R}^{m}$, we have

(11)
$\|v^{T}AR^{-1}\|_{2}(1+\varepsilon_{1})^{-1}\leq\|v^{T}A\hat{R}^{-1}\|_{2}\leq\|v^{T}AR^{-1}\|_{2}(1-\varepsilon_{1})^{-1}.$

The bound follows as $\kappa_{2}(AR^{-1})=\kappa_{2}(SQ_{A})$ (see Lemma 2.1).

Lemma 3.1 shows us, assuming $\varepsilon_{1}$ is sufficiently small, that
$\kappa_{2}(\hat{R})$ cannot be larger than about
$\kappa_{2}(SQ_{A})\kappa_{2}(A)$, and $\kappa_{2}(A\hat{R}^{-1})$ cannot be
much larger than $\kappa_{2}(SQ_{A})$.

As the next lemma shows, the error incurred by solving the triangular system
depends on the quantities in Lemma 3.1.

######  Lemma 3.3.

Let $\hat{R}^{-1}$ be the upper triangular matrix defined in Lemma 3.1 and let
$\hat{Y}$ be the matrix obtained by computing $A\hat{R}^{-1}$ with forward
substitution in finite precision. Define $\Delta Y=\hat{Y}-A\hat{R}^{-1}$ and
$\varepsilon_{2}=\|\Delta Y\|_{2}\|(A\hat{R}^{-1})^{\dagger}\|_{2}$. Assuming
that $\varepsilon_{2}<1$, we find that $\hat{Y}$ is full rank and

$\kappa_{2}(\hat{Y})\leq\frac{\kappa_{2}(A\hat{R}^{-1})+\varepsilon_{2}}{1-\varepsilon_{2}}.$

######  Proof 3.4.

This follows immediately from the relationship $\hat{Y}=A\hat{R}^{-1}+\Delta
Y$ and similar reasoning to the proof of Lemma 3.1.

In exact arithmetic, clearly $\kappa_{2}(Y)=\kappa_{2}(AR^{-1})$. The
numerical error incurred due to the triangular solve is controlled by
$\|\Delta Y\|$, which depends on $\kappa_{2}(\hat{R})$ and
$\kappa_{2}(A\hat{R}^{-1})$. We show that the error is small under the
assumption that $\kappa_{2}(A)u\ll 1$. As a result, the condition number of
$\hat{Y}$ will be bounded by a small constant related to the conditioning of
$S$ and $SQ_{A}$. Thus, the computed $AP$ has a reasonable condition number
provided that the numerical errors $\varepsilon_{1}$ and $\varepsilon_{2}$ are
small.

###  3.2 Bounding the numerical errors

We now bound the errors $\varepsilon_{1}$ and $\varepsilon_{2}$.

######  Lemma 3.5.

Assume the same setup as in Lemma 3.1 and Lemma 3.3, where
$\varepsilon_{1}=\|E\|_{2}\|R^{-1}\|_{2}$ and $\varepsilon_{2}=\|\Delta
Y\|_{2}\|(A\hat{R}^{-1})^{\dagger}\|_{2}$, and all quantities in finite
precision are computed with unit round-off $u$. Furthermore, define

(12)
$C_{1}:=\sqrt{sn}\gamma_{m}+\sqrt{n}\tilde{\gamma}_{sn}(1+\sqrt{sn}\gamma_{m}),\qquad
k(S):=\|S\|_{2}\|(SQ_{A})^{\dagger}\|_{2}.$

If $\kappa_{2}(A)<1/(C_{1}k(S))$, then we have

$\varepsilon_{1}\leq C_{1}k(S)\kappa_{2}(A),$

and

$\varepsilon_{2}\leq\frac{n\gamma_{n}\kappa_{2}(SQ_{A})(\kappa_{2}(A)\kappa_{2}(SQ_{A})+\varepsilon_{1})(1+\varepsilon_{1})}{\left(1-\varepsilon_{1}-\sqrt{n}\gamma_{n}(\kappa_{2}(A)\kappa_{2}(SQ_{A})+\varepsilon_{1})\right)(1-\varepsilon_{1})}.$

######  Proof 3.6.

We first note that $E=E_{1}+E_{2}$, where $E_{1}=\hat{B}-SA$ is the error
arising from matrix-matrix multiplication and $E_{2}=\tilde{Q}\hat{R}-\hat{B}$
is the error arising from the Householder QR algorithm. We can bound these
terms using Eq. 4 and Eq. 5 to find

$\displaystyle\|E\|_{2}$
$\displaystyle\leq\|E_{1}\|_{2}+\sqrt{n}\tilde{\gamma}_{sn}\|SA+E_{1}\|_{2}$
$\displaystyle\leq\sqrt{n}\tilde{\gamma}_{sn}\|S\|_{2}\|A\|_{2}+(1+\sqrt{n}\tilde{\gamma}_{sn})\sqrt{sn}\gamma_{m}\|S\|_{2}\|A\|_{2}=C_{1}\|S\|_{2}\|A\|_{2}.$

This is then combined with

$\sigma_{\min}(R)=\sigma_{\min}(QR)=\sigma_{\min}(SQ_{A}R_{A})\geq\sigma_{\min}(SQ_{A})\sigma_{\min}(A).$

As for $\|\Delta Y\|_{2}$, by Eq. 7 we have

$\|\Delta Y\|_{2}\leq\|\Delta
Y\|_{F}\leq\frac{\sqrt{n}\gamma_{n}\kappa_{2}(\hat{R})}{1-\sqrt{n}\gamma_{n}\kappa_{2}(\hat{R})}\|A\hat{R}^{-1}\|_{F}\leq\phi(\hat{R})\|A\hat{R}^{-1}\|_{2},$

where

(13)
$\phi(\hat{R}):=\frac{n\gamma_{n}\kappa_{2}(\hat{R})}{1-\sqrt{n}\gamma_{n}\kappa_{2}(\hat{R})}\leq\frac{n\gamma_{n}(\kappa_{2}(A)\kappa_{2}(SQ_{A})+\varepsilon_{1})}{1-\varepsilon_{1}-\sqrt{n}\gamma_{n}(\kappa_{2}(A)\kappa_{2}(SQ_{A})+\varepsilon_{1})}.$

We have $\varepsilon_{2}\leq\phi(\hat{R})\kappa_{2}(A\hat{R}^{-1})$, which can
be bounded using the result on $\kappa_{2}(A\hat{R}^{-1})$ in Lemma 3.1.

The condition that $\kappa_{2}(A)<1/(C_{1}k(S))$ is closely related to the
assumption $\kappa_{2}(A)\ll u^{-1}$ (see Section 3.3 for further discussion),
which ensures that $\varepsilon_{1}<\kappa_{2}(A)u\ll 1$. The bound
$\varepsilon_{2}\leq\phi(\hat{R})\kappa_{2}(A\hat{R}^{-1})$, where
$\phi(\hat{R})$ is defined in Eq. 13, tells us that the dominant term in the
bound for $\varepsilon_{2}$ is $n^{2}\kappa_{2}(SQ_{A})\kappa_{2}(A)u$, when
$\varepsilon_{1}<1$.

###  3.3 Bounding $\mathbf{\boldsymbol{\kappa}_{2}(\hat{Y})}$

We now combine Lemmas 3.1, 3.3 and 3.5 to obtain a bound on the condition
number of $\hat{Y}$, which is the computed version of $AR^{-1}$. We find that
$\kappa_{2}(\hat{Y})\leq 4\kappa_{2}(SQ_{A})+1$ provided that $\kappa_{2}(A)$
is sufficiently small.

######  Theorem 3.7.

Let $A\in\mathbb{R}^{m\times n}$ have linearly independent columns,
$S\in\mathbb{R}^{s\times n}$ have linearly independent rows, and $m>s>n$. Let
$A=Q_{A}R_{A}$ and $SA=QR$ be the economized QR factorizations of $A$ and
$B=SA$, respectively. Assume that

(14) $\kappa_{2}(SQ_{A})<49,\quad n^{2}u<1/201,\quad\frac{m}{n}+\sqrt{n}>200,$

and that $\kappa_{2}(A)$ is sufficiently small so that

(15) $\kappa_{2}(A)<\frac{1}{3C_{1}k(S)},$

where $C_{1}=\mathcal{O}((\sqrt{sn}m+sn^{3/2})u)$ and $k(S)$ are defined in
Eq. 12. Compute the preconditioned matrix $\hat{Y}$ as detailed in Lemma 3.1
and Lemma 3.3.333That is, compute the matrix-matrix product $SA$, compute its
$R$ factor by the Householder QR algorithm, and finally compute the
preconditioned matrix by solving $AR^{-1}$ with forward substitution. Then,

$\kappa_{2}(\hat{Y})\leq 4\kappa_{2}(SQ_{A})+1.$

######  Proof 3.8.

Throughout the proof, we assume the notation introduced in Lemmas 3.1, 3.3 and
3.5. By assumption Eq. 15 we immediately see that $\varepsilon_{1}<1/3$,
defined in Lemma 3.1, and hence $\hat{R}$ and $A\hat{R}^{-1}$ are full rank.
It follows that $\kappa_{2}(A\hat{R}^{-1})\leq 2\kappa_{2}(SQ_{A})$.
Furthermore, also by Lemma 3.1, we find

(16)
$\kappa_{2}(\hat{R})\leq\frac{\kappa_{2}(R)+\varepsilon_{1}}{1-\varepsilon_{1}}\leq\frac{1}{2}(3\kappa_{2}(SQ_{A})\kappa_{2}(A)+1)<\frac{1}{2}\left(\frac{1}{C_{1}}+1\right).$

By using the fact that $m>s>n$, we can show that
$C_{1}>(m+n^{3/2})\gamma_{n},$ which can be substituted into Eq. 16. One can
use the bound on $\kappa_{2}(\hat{R})$ to then bound $\phi(\hat{R})$, defined
in Eq. 13, as

$\phi(\hat{R})=\frac{n\gamma_{n}\kappa_{2}(\hat{R})}{1-\sqrt{n}\gamma_{n}\kappa_{2}(\hat{R})}<\frac{\frac{1}{2}n\left(\frac{1}{m+n^{3/2}}+\gamma_{n}\right)}{1-\frac{1}{2}\sqrt{n}\left(\frac{1}{m+n^{3/2}}+\gamma_{n}\right)}=\frac{\sqrt{n}c(m,n)}{1-c(m,n)},$

where

$c(m,n)=\frac{1}{2}\sqrt{n}\left(\frac{1}{m+n^{3/2}}+\gamma_{n}\right).$

The assumptions in Eq. 14 imply that
$c(m,n)<(4\sqrt{n}\kappa_{2}(SQ_{A})+1)^{-1}$, where we used that
$n^{2}u<1/201$, i.e. $n\gamma_{n}<1/200$. Now we have that
$\phi(\hat{R})<(4\kappa_{2}(SQ_{A}))^{-1}$ and so

$\varepsilon_{2}\leq\phi(\hat{R})\kappa_{2}(A\hat{R}^{-1})\leq\kappa_{2}(SQ_{A})\phi(\hat{R})\frac{1+\varepsilon_{1}}{1-\varepsilon_{1}}<\frac{1}{2}.$

By Lemma 3.3, $\hat{Y}$ is full rank and its condition number is bounded by

$\kappa_{2}(\hat{Y})\leq\frac{\kappa_{2}(A\hat{R}^{-1})+\varepsilon_{2}}{1-\varepsilon_{2}}<4\kappa_{2}(SQ_{A})+1.$

It is worth discussing the assumptions of Theorem 3.7. Consider the assumption
in Eq. 15. The $k(S)$ term (see Eq. 12) is a small constant depending on the
particulars of the embedding matrix. The constant $C_{1}$ (see Eq. 12) can be
bounded by the product of a low-degree polynomial in $m$, $n$, and $s$ and the
unit-round-off $u$, say $C_{1}\leq p(m,n,s)u$. We specifically have
$C_{1}=\mathcal{O}(\sqrt{sn}(m+\sqrt{s}n)u)$, as $\gamma_{m}=\mathcal{O}(mu)$
and $\tilde{\gamma}_{sn}=\mathcal{O}(snu)$. In classical stability analysis,
this polynomial term $p(m,n,s)$ can be more or less ignored. As a result, a
condition of the form $\kappa_{2}(A)<(p(m,n,s)u)^{-1}$ is often loosely
restated as $\kappa_{2}(A)\ll u^{-1}$.

The rationale behind ignoring these $m$, $n$, and $s$ terms is that classic
stability analysis, as we performed in this section, provides us with worst-
case bounds on numerical errors. These are generally pessimistic [11];
composing $n$ operations would only result in an error of order $nu$ if each
rounding error is of the same sign and of maximum magnitude [12]. We could
improve these bounds by using probabilistic backward error analysis, which
would give us results proportional to $\sqrt{n}u$ instead of $nu$ [12, 13].
However, it must be noted that even probabilistic stability analysis would
theoretically result in a large factor multiplied by $u$.

The other assumptions of Theorem 3.7, in Eq. 14, are reasonable in practice.
The condition number of $SQ_{A}$ is generally observed to be bounded by a
modest number, such as $5$ or $10$, with high probability. Furthermore, we
assume that the problem dimensions are small compared to the unit round-off.
Most importantly, Theorem 3.7 informs us that, provided that $\kappa_{2}(A)\ll
u^{-1}$, the condition number of the preconditioned matrix is independent of
$\kappa_{2}(A)$ (with high probability).

###  3.4 The backward error of sketch-and-apply Blendenpik

Until now, we have only considered the condition number of the preconditioned
matrix $AP$. The analysis in this section related rounding errors incurred in
the final steps — solving $(AR^{-1})z=b$ with LSQR and computing $x=R^{-1}z$
with back substitution — to the backward error of sketch-and-apply. It is not
immediate that the last step in the algorithm preserves the numerical
stability of the algorithm. After all, solving linear systems with $R$ in LSQR
iterations is one potential reason for the numerical instability of sketch-
and-precondition Blendenpik. Nonetheless, we show that sketch-and-apply
Blendenpik performs as well as LSQR on very well-conditioned matrices and that
the numerical error in the last step is $\mathcal{O}(u)$ instead of
$\mathcal{O}(u\kappa_{2}(A))$ (see Theorem 3.9).

We show that the backward error in sketch-and-apply Blendenpik is of the same
order as the backward error from unpreconditioned LSQR with
$\text{fl}(AR^{-1})$, under the assumptions of Theorem 3.7.

######  Theorem 3.9.

Assume the embedding matrix $S\in\mathbb{R}^{s\times m}$ and
$A\in\mathbb{R}^{m\times n}$ satisfy assumptions Eq. 14 and Eq. 15 in Theorem
3.7. Apply Algorithm 2 and assume that LSQR terminates at an iterate $\hat{z}$
that satisfies the backward error444Note that the backward error term
$\Delta\hat{Y}$ is not the same as the numerical error $\Delta Y$ resulting
from solving the triangular system $Y=A\hat{R}^{-1}$ in finite precision.

(17) $\hat{z}=(\hat{Y}+\Delta\hat{Y})^{\dagger}(b+\delta b).$

Then, the computed solution $\hat{x}$ to $Ax=b$ has the backward error
$(A+\Delta A)\hat{x}=b+\delta b$ satisfying

$\|\Delta
A\|_{2}<\|S\|_{2}\|A\|_{2}\left(6.04n\gamma_{n}\|(SQ_{A})^{\dagger}\|_{2}+2.01\|\Delta\hat{Y}\|_{2}\right).$

######  Proof 3.10.

The final step of Algorithm 2 involves solving the triangular system
$\hat{R}x=\hat{z}$. We know from Eq. 8 that the computed solution $\hat{x}$
satisfies

$(\hat{R}+\Delta\hat{R})\hat{x}=\hat{z},\quad\|\Delta\hat{R}\|_{2}\leq\gamma_{n}\sqrt{n}\|\hat{R}\|_{2}.$

We now have

$\hat{x}=(\hat{R}+\Delta\hat{R})^{-1}\hat{z}=(\hat{R}+\Delta\hat{R})^{-1}(\hat{Y}+\Delta\hat{Y})^{\dagger}(b+\delta
b)=(A+\Delta A)^{\dagger}(b+\delta b),$

where $\Delta
A=(\hat{Y}\hat{R}-A)+\Delta\hat{Y}\hat{R}+\hat{Y}\Delta\hat{R}+\Delta\hat{Y}\Delta\hat{R}.$
For the $\hat{Y}\hat{R}-A$ term we have that each row satisfies

$\hat{y}_{i}^{T}\hat{R}-a_{i}^{T}=-\hat{y}_{i}^{T}\Delta_{i}\hat{R},\qquad|\Delta_{i}\hat{R}|\leq\gamma_{n}|\hat{R}|,$

so that $\|\hat{Y}\hat{R}-A\|_{2}\leq
n\gamma_{n}\|\hat{R}\|_{2}\|\hat{Y}\|_{2}$. It follows that

$\|\Delta
A\|_{2}\leq\|\hat{R}\|_{2}\left((n+\sqrt{n})\gamma_{n}\|\hat{Y}\|_{2}+(1+\sqrt{n}\gamma_{n})\|\Delta\hat{Y}\|_{2}\right),$

where $\|\hat{R}\|_{2}\leq(1+C_{1})\|S\|_{2}\|A\|_{2}<2\|S\|_{2}\|A\|_{2},$
and

$\|\hat{Y}\|_{2}\leq\|A\hat{R}^{-1}\|_{2}+\|\Delta
Y\|_{2}\leq\left(1+\phi(\hat{R})\right)\frac{\|(SQ_{A})^{\dagger}\|_{2}}{1-\varepsilon_{1}}<1.51\|(SQ_{A})^{\dagger}\|_{2},$

where the last inequality follows from the assumptions Eq. 14 and Eq. 15. We
obtain

$\|\Delta A\|\leq
2\|S\|_{2}\|A\|_{2}\left(3.02n\gamma_{n}\|(SQ_{A})^{\dagger}\|_{2}+1.005\|\Delta\hat{Y}\|_{2}\right),$

where we used that $n+\sqrt{n}\leq 2n$ and $\sqrt{n}\gamma_{n}\leq
n\gamma_{n}\leq 1/200$ by Eq. 14.

The moral of Theorem 3.9 is that sketch-and-apply Blendenpik is backward
stable as long as LSQR on $\hat{Y}$ is backward stable, that is,
$\|\Delta\hat{Y}\|_{2}$ and $\|\delta b\|_{2}/\|b\|_{2}$ are both $O(u)$ in
Eq. 17. Under the conditions of Theorem 3.9, the matrix $\hat{Y}$ is well-
conditioned with high probability. The backward stability of LSQR on a well-
conditioned matrix is not an undisputedly clear fact. Practical evidence is
abundant and several studies [3, 7, 8, 9, 17] have addressed the convergence
of the conjugate gradient (CG) method (LSQR is a stable implementation of CG
applied to the normal equation), we are unaware of a precise result that
proves CG (or LSQR) applied to a well-conditioned system is backward stable.
Such stability analysis would depend on the precise implementation of the
algorithm (see Section 6).

Remarkably, the final step of the algorithm does not influence the accuracy of
the solution, i.e., solving $Rx=y$ with back substitution. This is despite the
fact that $R$ is ill-conditioned when $A$ is. The mechanism with which
stability is established is similar to the backward stability of an algorithm
based on repeated CholeskyQR [27].

##  4 On the numerical stability of smoothed sketch-and-apply Blendenpik

We now analyze smoothed sketch-and-apply Blendenpik (see Algorithm 3). When
$\kappa_{2}(A)\ll 1/u$, Theorem 3.7 states that sketch-and-apply computes a
preconditioned matrix with a condition number independent of $\kappa_{2}(A)$
with high probability. Moreover, we can reliably estimate $\kappa_{2}(A)$ by
$\kappa_{2}(SA)$. We do not have a guarantee on the stability of sketch-and-
apply Blendenpik if $\kappa_{2}(A)\gtrsim 1/u$. It is in this context that we
consider smoothing.

One could consider Algorithm 3 from two points of view. Firstly, a
practitioner may add an $\mathcal{O}(u\|A\|_{2})$ Gaussian perturbation to a
severely ill-conditioned problem. As a result, the perturbed problem can be
accurately solved with sketch-and-apply. Secondly, floating point errors
arising from the representation of any matrix in finite precision could be
modeled as a Gaussian perturbation and then smoothing partially explains the
success of sketch-and-apply to highly ill-conditioned LS problems.

###  4.1 Bounding $\mathbf{\boldsymbol{\kappa}_{2}(\hat{Y})}$

Similar to the stability analysis in Section 3, we have the following steps in
finite precision after we have smoothed $A$:

  1. 1.

$\tilde{A}=A+\sigma G/\sqrt{m}$,

  2. 2.

$\hat{B}=S\tilde{A}+E_{1}$,

  3. 3.

$\tilde{Q}\hat{R}=\hat{B}+E_{2}$,

  4. 4.

$\hat{Y}=\tilde{A}\hat{R}^{-1}+\Delta Y$,

where the hats denote that we are accounting for rounding errors and the
scaling parameter $\sigma$ controls how much $A$ is perturbed. Assuming
$\sigma$ is sufficiently large to ensure that $\tilde{A}$ satisfies Eq. 15, we
then show that the stability analysis from Section 3 carries over to
$\tilde{A}$.

The following theorem tells us how large $\sigma$ needs to be.

######  Theorem 4.1.

Let $A\in\mathbb{R}^{m\times n}$ with $m\geq 3n$, $S\in\mathbb{R}^{s\times n}$
have linearly independent rows, and let $m>s>n$. Set

(18) $\tilde{A}=A+\sigma G/\sqrt{m},\qquad\sigma=52\|A\|_{2}k(S)sm\sqrt{n}u,$

where $G\in\mathbb{R}^{m\times n}$ is a standard Gaussian matrix with i.i.d.
entries, $k(S)$ is defined in Eq. 12 and $u$ is the unit round-off. Let
$\tilde{A}=Q_{\tilde{A}}R_{\tilde{A}}$ and $S\tilde{A}=QR$ be economized QR
factorizations of $\tilde{A}$ and $B=S\tilde{A}$, respectively. Lastly, assume
that555We also assume that $cn<m$ for the same integer $c$ defined in
$\tilde{\gamma}_{sn}$ to simplify the notation.

(19) $\kappa_{2}(SQ_{\tilde{A}})<49,\qquad\frac{m}{n}+\sqrt{n}>200,\qquad
smu<1/67.$

Now compute first the product $S\tilde{A}$, then the upper triangular factor
$R$ by the Householder QR algorithm on $S\tilde{A}$, and finally compute
$\hat{Y}=\tilde{A}R^{-1}$ with forward substitution. Assume all computations
are performed in finite precision. Then,

$\mathbb{P}\left(\kappa_{2}(\hat{Y})>4\kappa_{2}(SQ_{\tilde{A}})+1\right)<2^{n-m}.$

######  Proof 4.2.

We start by showing that $\sigma$ in Eq. 18 is sufficiently large to ensure
that $\tilde{A}$ satisfies Eq. 15 in Theorem 3.7. To this end, note the
assumptions in Eq. 19 imply that $mu<1/64$ and $csnu<1/64$ so that we have
$\gamma_{m}<1.02mu$ and $\tilde{\gamma}_{sn}<1.02csnu<1.02smu$. We can then
bound $C_{1}$ defined in Eq. 12 as

$C_{1}<1.02\times 65sm(1+\sqrt{n})u/64<2.1sm\sqrt{n}u.$

As a result

$\mathbb{P}\left(\kappa_{2}(\tilde{A})>\frac{1}{3C_{1}k(S)}\right)<\mathbb{P}\left(\kappa_{2}(\tilde{A})>\frac{8.25}{52k(S)sm\sqrt{n}u}\right)<2^{n-m}.$

The final result on $\kappa_{2}(\hat{Y})$ follows by applying Theorem 3.7 to
$\tilde{A}$, which is possible since the last assumption in Eq. 19 is stronger
than $nu<1/201$.

Theorem 4.1 shows that the magnitude of the perturbation $\sigma G/\sqrt{m}$
needs to be approximately $\mathcal{O}(sm\sqrt{n}u)$ to ensure the condition
number of $\tilde{A}$ is sufficiently small, that is, $\kappa_{2}(A)<1/u$. The
low-degree polynomial term $sm\sqrt{n}$ arises from Eq. 15, which includes the
constant $C_{1}=\mathcal{O}(\sqrt{sn}(m+\sqrt{s}n)u)$. As discussed in Section
3.4, these are often pessimistic bounds. Probabilistic error analysis provides
a theoretical justification to select $\sigma$ proportional to
$\sqrt{sm}n^{1/4}u$. This may still be a relatively large perturbation if the
dimensions of the problem are enormous, enlarging the backward error to an
unacceptable level (see Section 4.2).

###  4.2 The backward error of smoothed sketch-and-apply Blendenpik

In Section 4.1, we showed that a small additive random perturbation to a
small-skinny matrix ensures that the condition number of the perturbed matrix
is sufficiently small. This allows us to conclude that $\kappa_{2}(\hat{Y})$
is small with high probability. We now show that smoothed sketch-and-apply
Blendenpik (see Algorithm 3) computes a backward error similar to Theorem 3.9
with an additional term from the additive perturbation.

######  Theorem 4.3.

Assume the same setup as in Theorem 4.1. Apply Algorithm 2 to $\tilde{A}x=b$
and assume LSQR terminates at an iterate $\hat{z}$ that satisfies the backward
error

$\hat{z}=(\hat{Y}+\Delta\hat{Y})^{\dagger}(b+\delta b).$

Then, the computed solution $\hat{x}$ to $Ax=b$ has the backward error
$(A+\Delta A)\hat{x}=b+\delta b$ satisfying

(20) $\|\Delta
A\|_{2}\leq\|A\|_{2}\left[\vphantom{\frac{1}{2}}\|S\|_{2}\left(1+197k(S)sm\sqrt{n}u\right)\right.\\\
\left.\left(6.04n\gamma_{n}\|(SQ_{\tilde{A}})^{\dagger}\|_{2}+2.01\|\Delta\hat{Y}\|_{2}\right)+197k(S)sm\sqrt{n}u\vphantom{\frac{1}{2}}\right]$

with probability at least $1-2^{n-m}-10^{-3}$.

######  Proof 4.4.

We combine Theorem 3.9 and Theorem 4.1 to obtain the result. Additionally, we
use the following classic result in the analysis of Gaussian matrices [6]:

$\mathbb{P}\left(\|G\|_{2}/\sqrt{m}>1+\sqrt{\frac{n}{m}}+\frac{t}{\sqrt{m}}\right)\leq
e^{-t^{2}/2},$

applied with $t=3.8$. This shows
$\|\tilde{A}\|_{2}\leq\|A\|_{2}+\sigma\|G\|_{2}/\sqrt{m}\leq\|A\|_{2}+3.78\sigma$
with probability at least $10^{-3}$, where
$\sigma=52\|A\|_{2}k(S)sm\sqrt{n}u$.

Theorem 4.3 allows us to conclude that under mild assumptions, smoothed
sketch-and-apply Blendenpik results in a backward error of the same order as
the backward error of unpreconditioned LSQR on $\hat{Y}z=b$. Under the
conditions of the theorem, $\hat{Y}$ is indeed a well-conditioned matrix and
LSQR is expected to obtain backward stable solutions. We removed the
assumptions on $\kappa_{2}(A)$.

Theorem 4.3 suggests that the perturbation needs to be of size
$\mathcal{O}(sm\sqrt{n}u)$. This would be impermissible for many problems,
especially if the dimensions are large. We recommend a perturbation of
magnitude approximately $10\|A\|_{2}u$. Although we lack the theoretical
evidence for this choice, it works in our experiments (see Section 5.2).

Smoothing a LS problem should be done with much caution. There are many
contexts where $\kappa_{2}(A)>u^{-1}$, yet perturbing the problem is
unnecessary. For instance, if the ill-conditioning is caused by poor column
scaling, sketch-and-apply will still compute accurate solutions without
smoothing (see Section 6). In this sense, $\kappa_{2}(A)$ is not an effective
predictor of whether smoothing will be beneficial. We suggest one should
smooth only if sketch-and-apply does not converge (see Algorithm 4).

##  5 Master algorithm and numerical experiments

We now present the implementation of the sketch-and-apply algorithms (see
Algorithm 4) in a similar fashion to the overall solver in Blendenpik [1]. We
perform the standard initial steps, i.e., drawing an embedding matrix $S$,
computing the sketch $SA$, and computing the upper triangular factor $R$ of
the QR factorization of $SA$ using the Householder QR algorithm. Next, we
compute the preconditioned matrix $Y=AR^{-1}$ explicitly with forward
substitution and use (unpreconditioned) LSQR to solve $Yz=b$. If LSQR
converges to the desired tolerance, we compute $x=R^{-1}z$ with back
substitution. In this case, we assume our algorithm computed a backward stable
solution.

In the case where LSQR does not converge, we perturb/smooth $A$. We compute
$\tilde{A}=A+\sigma G/\sqrt{m}$ for a Gaussian matrix $G$ with
$\sigma=10\|A\|_{2}u$ and perform the sketch-and-apply steps on
$\tilde{A}x=b$. We recommend doing this because there are examples, although
rare, where sketch-and-apply Blendenpik converges to an accurate solution only
after smoothing (see Section 5.2). The spectral norm of $A$ can be estimated,
for instance, with the sketch $SA$ or using the power method.

Algorithm 4 Master algorithm for sketch-and-apply to solve Eq. 1. Here, HHQR
refers to the Householder QR algorithm.

1: Draw a random sketching matrix $S\in\mathbb{R}^{s\times m}$, where $m\gg
s>n$

2: Compute $B=SA$

3: Compute the triangular factor $R$ of a QR factorization of $B$ using HHQR

4: Compute $Y=AR^{-1}$ with forward substitution

5: Solve $Yz=b$ with LSQR and no preconditioner

6: if LSQR converges to the desired tolerance then

7: Compute $x=R^{-1}z$ with back substitution

8: else

9: Draw a random standard Gaussian matrix $G\in\mathbb{R}^{m\times n}$ with
i.i.d. entries

10: Compute $\tilde{A}=A+\sigma G/\sqrt{m}$ for $\sigma=10\|A\|_{2}u$

11: Compute $B=S\tilde{A}$

12: Compute the triangular factor $R$ of a QR factorization of $B$ using HHQR

13: Compute $Y=\tilde{A}R^{-1}$ with forward substitution

14: Solve $Yz=b$ with LSQR and no preconditioner

15: Compute $x=R^{-1}z$ with back substitution

16: end if

It should be noted that without smoothing, LSQR iteration can converge to a
solution with a sub-optimal residual (see Fig. 6 (left)). However, it is
infeasible to compute the backward error when $m$ is large as it requires
computing the smallest singular value of an $m\times(m+n)$ matrix. If the
residual obtained from sketch-and-apply is less accurate than desired, one
could estimate the condition number of $A$ via the condition number of $R$. If
the condition number estimate of $R$ is of the same order of magnitude as
$u^{-1}$, then the problem is severely ill-conditioned and one could try
smoothing to potentially improve the accuracy of the computed solution.
Smoothing can improve the convergence rate of LSQR, even if sketch-and-apply
obtains a backward stable solution without smoothing (see Fig. 5). Nonetheless
one should be careful with smoothing because examples exist where
$\kappa_{2}(A)$ is large but smoothing remains unnecessary (see Fig. 5).
Exactly when and when not to smooth remains unclear.

###  5.1 Numerical experiments with sketch-and-apply Blendenpik

Now we examine the performance of sketch-and-precondition Blendenpik, sketch-
and-apply Blendenpik, and the QR-based direct solver. The first experiment
considers ill-conditioned LS problems that are randomly generated (see Section
1.1 for details on how $A$ and $b$ are formed) and can be stored in local
cache. The execution times are shown in Fig. 3. We find that on this smaller
scale, sketch-and-apply Blendenpik is competitive relative to QR in terms of
computational time. While sketch-and-precondition is the computationally most
efficient, it fails to compute a backward stable solution (see Fig. 1).

Figure 3: The execution time to solve large LS problems with sketch-and-
precondition Blendenpik, sketch-and-apply Blendenpik, and Householder QR. For
each problem $\kappa(A)=10^{10}$. Left: Timings for $n=4000$ and $2^{12}\leq
m\leq 2^{20}$. Right: Timings for $m=5\times 10^{5}$ and $2^{7}\leq n\leq
2^{12}$.

Next, we investigate the computational performance on LS problems for which
the $m\times n$ matrix $A$ is not stored in local cache (see Fig. 4). Instead,
it is stored on disk in a CSV format. Working with $A$ in this form requires
the algorithms to be implemented again. For example, the QR algorithm now
reads a tall-skinny matrix ten columns at a time, performs Householder
reflectors column-by-column, before writing the vector representing the
Householder reflector to disk. The column is deleted from memory and the next
column is requested from disk. In this way, the $n\times n$ matrix $R$ is
constructed column-by-column and only one pass over the matrix $A$ is
required, although it is necessary to pass over the CSV file that stores the
Householder reflectors many times. The implementation for sketch-and-
precondition is easy to modify as one only needs to implement matrix-vector
products with $A$ and $A^{\top}$ when $A$ is on disk. For this, we read ten
columns of $A$ (or ten rows of $A^{T}$) into local cache before partially
computing the matrix-vector product, deleting the columns, and then requesting
the next ten cols (or rows). Sketch-and-apply is the same except we need to
compute the preconditioned matrix $\hat{Y}$, which we write to a separate CSV
file on disk. The Householder QR algorithm is computationally more efficient
for small dimensions than sketch-and-precondition and sketch-and-apply.
However, for large LS problems, sketch-and-precondition and sketch-and-apply
start to become more computationally efficient. In Fig. 4, LSQR in sketch-and-
apply requires about 20 iterations, independent of $n$. Hence, the number of
times that $A$ has to be streamed from disk to local cache is modest.

Figure 4: The execution time to solve large-scale LS problems of varying sizes
with sketch-and-precondition Blendenpik, sketch-and-apply Blendenpik, and
Householder QR. These matrices are stored in a CSV file on disk rather than in
local cache memory. For each problem $\kappa_{2}(A)=10^{15}$, $m=10^{4}$, and
$10\leq n\leq 10^{3}$. We set the noise level to zero and the tolerance in
LSQR to $10^{-10}$.

###  5.2 Numerical experiments with smoothed sketch-and-apply Blendenpik

Sketch-and-apply Blendenpik, its smoothed version, and Householder QR can
compute accurate solutions to LS problems, even when $\kappa_{2}(A)>1/u$ (see
Fig. 5, the problems are generated as explained in Section 1.1). Therefore,
smoothing is only sometimes required for extremely ill-conditioned least
squares problems. However, even when smoothing is unnecessary, there is some
potential benefit to smoothing because the perturbed LS problem can lead to
rapid LSQR convergence. Of course, there is a trade-off here as one converges
rapidly to an accurate solution of the perturbed problem, and the computed
solution may be less accurate for the original LS problem of interest.

Figure 5: Smoothing is not always needed for LS problems with
$\kappa_{2}(A)>1/u$. These matrices are formed by drawing singular vectors
from the Haar distribution and letting the singular values decay exponentially
from 1 to $1/\kappa_{2}(A)$. For comparison, we show the residual error
computed by Householder QR.

However, there are also examples for which smoothing provides a more accurate
solution to the original problem of $\min_{x}\|Ax-b\|$, not
$\min_{x}\|\tilde{A}x-b\|$ (see Fig. 6). For example, we take an LS problem
involving a $1000\times 100$ Kahan matrix with $\theta=1.1$ from the MATLAB
gallery collection and another one involving a column-scaled $1000\times 10$
Vandermonde matrix involving equally spaced points between $-1$ and $1$. The
Kahan and Vandermonde matrices are designed so that their condition numbers
are $>1/u$. Without smoothing, even the QR-based algorithm fails to compute an
accurate solution along with all the other methods. This is because the
problem is so ill-conditioned that even a backward stable solution can behave
wildly. However, after smoothing, sketch-and-apply Blendenpik computes an
accurate solution to the original LS problems. It should be noted that one
could also smooth before solving the problems with QR to obtain more accurate
solutions.

Figure 6: Smoothing can solve extremely ill-conditioned LS problems, even when
sketch-and-apply Blendenpik and Householder QR cannot. For the Vandermonde
matrix, the residuals when using QR or sketch-and-apply Blendenpik both blow-
up. Specifically, the relative residual with the QR-based solver is about
$10^{115}$, and it is $10^{85}$ for sketch-and-apply Blendenpik.

##  6 Discussion

We have shown that sketch-and-precondition algorithms, such as Blendenpik [1],
are numerically unstable in their standard form for solving LS problems. We
have stabilized the algorithm by explicitly computing the preconditioned
matrix $AP$ and using an unpreconditioned iterative solver on $AP$. We coined
this modification sketch-and-apply.

###  6.1 The numerical stability of data-driven preconditioners

The numerical instabilities observed in sketch-and-precondition Blendenpik
raise larger questions on the stability of iterative methods with data-driven
preconditioners for LS problems. Here, we refer to a data-driven
preconditioner as a preconditioner constructed directly from $A$, without
knowing where $A$ came from (such as the discretization of a continuous
problem). Fig. 2 shows that even using $R_{A}$ (where $A=Q_{A}R_{A}$) as a
preconditioner—the perfect data-driven preconditioner—does not lead to a
backward stable solution. We suspect this is due to the numerical errors
incurred each time the ill-conditioned preconditioner is applied in an
iterative solver. Is it possible for a data-driven preconditioner to avoid
compounding these rounding errors?

###  6.2 The numerical stability of iterative least squares solvers

As to the numerical stability of sketch-and-apply, we were not able to state
that the backward error $\|\Delta A\|_{2}$ is $\mathcal{O}(u\|A\|_{2})$ (see
Theorems 3.9 and 4.3). Instead, we have shown that $\|\Delta
A\|_{2}/\|A\|_{2}$ is of the same order as the backward error
$\|\Delta\hat{Y}\|$ incurred when $\hat{Y}z=b$ is solved with unpreconditioned
LSQR, where $\hat{Y}$ is well-conditioned. It has proven challenging to
understand the literature on the numerical stability of CG-like iterative
solvers such as LSQR. Various works seem to strongly hint at backward
stability under assumptions on the condition number, but use computational
results to complete the claim, see [3, 8]. The numerical stability depends
strongly on the specific implementation used in a way that we are yet to
understand fully.

###  6.3 Variants of sketch-and-precondition

We have assumed specific choices for how the randomized preconditioner is
constructed. In LSRN [16], for instance, the preconditioner is chosen to be
$P=V\Sigma^{-1}$, where $SA=U\Sigma V^{T}$ is the SVD of the sketch. Numerical
experiments show us that SVD-based sketch-and-precondition methods incur
similar numerical instabilities as the QR-based variant. The sketch-and-apply
technique can also be used with LSRN ideas where the preconditioned matrix is
computed as $Y=AV\Sigma^{-1}$. We suspect one can also prove similar results
to our sketch-and-apply Blendenpik analysis for sketch-and-apply LSRN.
Versions of Lemmas 3.1 and 3.3 hold almost identically for LSRN, with $R$
replaced by $\Sigma V$. The key difference is the numerical errors
$\varepsilon_{1}$ and $\varepsilon_{2}$, which we have not investigated. Since
computing the SVD is typically more expensive than a QR factorization, we
recommend sketch-and-apply Blendenpik. One notable exception is when solving a
sequence of Tikhonov regularized LS problems for various regularization
parameters [15, 16].

###  6.4 Ill-conditioning caused by poor column scaling

If the ill-conditioning in an LS problem is caused by poor column scaling in
$A$, we strongly recommend using QR-based sketch-and-apply techniques. The
reason is that computing $AR^{-1}$ is invariant to column scaling, while the
SVD is not. In fact, the assumption that $\kappa_{2}(A)\ll u^{-1}$ in Section
3 can be replaced by

$\min_{D}\\{\kappa_{2}(AD)\,:\,D=\text{diag($b^{k_{i}}$),
$k_{i}\in\mathbb{Z},\,i=1,\dots,n$}\\}\ll u^{-1},$

where $b$ is the machine base (usually $b=2$). Of course, one can also pre-
process an LS problem by scaling the columns of $A$ by powers of the machine
base so that each column has a norm that is close to $1$.

## Acknowledgments

We thank Erin Carson and Zdeněk Strakoš for their valuable comments regarding
the numerical stability of LSQR. We thank Françoise Tisseur for her input on
Blendenpik in finite precision arithmetic. We are indebted to Ilse Ipsen and
Michael Mahoney for their presentation and subsequent discussions on LS
problems, which occurred during the “Complexity of Matrix Computations”
seminar on 1st September 2021.

## References

  * [1] H. Avron, P. Maymounkov, and S. Toledo.  Blendenpik: Supercharging LAPACK’s least-squares solver.  SIAM J. Sci. Comput., 32(3):1217–1236, 1 2010.
  * [2] Å. Björck.  Numerical Methods for Least Squares Problems.  SIAM, 1996.
  * [3] Å. Björck, T. Elfving, and Z. Strakoš.  Stability of conjugate gradient and Lanczos methods for linear least squares problems.  SIAM J. Mat. Anal. Appl., 19(3):720–736, 1998.
  * [4] C. Boutsidis and P. Drineas.  Random projections for the nonnegative least-squares problem.  Lin. Alg. Appl., 431(5-7):760–771, 2009.
  * [5] P. Bürgisser and F. Cucker.  Smoothed analysis of Moore–Penrose inversion.  SIAM J. Mat. Anal. Appl., 31(5):2769–2783, 2010.
  * [6] K. R. Davidson and S. J. Szarek.  Local operator theory, random matrices and Banach spaces.  Handbook of the geometry of Banach spaces, 1(317-366):131, 2001\.
  * [7] A. Greenbaum.  Behavior of slightly perturbed Lanczos and conjugate-gradient recurrences.  Lin. Alg. Appl., 113:7–63, 1989.
  * [8] A. Greenbaum.  Estimating the attainable accuracy of recursively computed residual methods.  SIAM J. Mat. Anal. Appl., 18(3):535–551, 1997.
  * [9] A. Greenbaum and Z. Strakos.  Predicting the behavior of finite precision Lanczos and conjugate gradient computations.  SIAM J. Mat. Anal. Appl., 13(1):121–137, 1992.
  * [10] N. Halko, P.-G. Martinsson, and J.A. Tropp.  Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions.  SIAM Rev., 53(2):217–288, 2011.
  * [11] N. J. Higham.  Accuracy and Stability of Numerical Algorithms.  SIAM, 2002.
  * [12] N. J. Higham and T. Mary.  A new approach to probabilistic rounding error analysis.  SIAM J. Sci. Comput., 41(5):A2815–A2835, 2019.
  * [13] N. J. Higham and T. Mary.  Sharper probabilistic backward error analysis for basic linear algebra kernels with random data.  SIAM J. Sci. Comput., 42(5):A3427–A3446, 2020.
  * [14] P.-G. Martinsson and J. A. Tropp.  Randomized numerical linear algebra: Foundations and algorithms.  Acta Numer., 29:403–572, 2020.
  * [15] M. Meier and Y. Nakatsukasa.  Randomized algorithms for Tikhonov regularization in linear least squares.  arXiv:2203.07329, 2022.
  * [16] X. Meng, M. A. Saunders, and M. W. Mahoney.  LSRN: A parallel iterative solver for strongly over- or underdetermined systems.  SIAM J. Sci. Comput., 36(2), 2014.
  * [17] G. Meurant and Z. Strakoš.  The Lanczos and conjugate gradient algorithms in finite precision arithmetic.  Acta Numer., 15:471–542, 2006.
  * [18] C. C. Paige and M. A. Saunders.  LSQR: An algorithm for sparse linear equations and sparse least squares.  ACM Trans. Math. Soft., 8(1):43–71, 1982.
  * [19] V. Rokhlin and M. Tygert.  A fast randomized algorithm for overdetermined linear least-squares regression.  Proc. Nat. Acad. Sci., 105(36):13212–13217, 2008.
  * [20] J. A. Tropp.  Improved analysis of the subsampled randomized Hadamard transform.  Adv. Adapt. Data Anal., 3(01n02):115–126, 2011.
  * [21] J. A. Tropp.  Randomized block Krylov methods for approximating extreme eigenvalues.  Numer. Math., 150(1):217–255, 2022.
  * [22] R. Vershynin.  Introduction to the non-asymptotic analysis of random matrices.  arXiv:1011.3027, 2010.
  * [23] B. Waldén, R. Karlson, and J.-G. Sun.  Optimal backward perturbation bounds for the linear least squares problem.  Numer. Lin. Alg. Appl., 2(3):271–286, 1995.
  * [24] A. J. Wathen and T. Rees.  Chebyshev semi-iteration in preconditioning for problems including the mass matrix.  Electron. Trans. Numer. Anal, 34:125–135, 2008.
  * [25] P.-Å. Wedin.  Perturbation theory for pseudo-inverses.  BIT Num. Math., 13(2):217–232, 1973.
  * [26] D. P. Woodruff.  Sketching as a tool for numerical linear algebra.  Foundations and Trends® in Theoretical Computer Science, 10(1–2):1–157, 2014.
  * [27] Y. Yamamoto, Y. Nakatsukasa, Y. Yanagisawa, and T. Fukaya.  Roundoff error analysis of the Cholesky QR2 algorithm.  Electron. Trans. Numer. Anal, 44:306–326, 2015.