A sketch-and-select Arnoldi processS. Güttel and I. Simunec

# A sketch-and-select Arnoldi process

Stefan Güttel Department of Mathematics, The University of Manchester, M13 9PL
Manchester, United Kingdom<EMAIL_ADDRESS>S. G.
acknowledges a fellowship from The Alan Turing Institute under the EPSRC grant
EP/W001381/1 and a Royal Society Industry Fellowship IF/R1/231032. Igor
Simunec Scuola Normale Superiore, 56126 Pisa, Italy<EMAIL_ADDRESS>
###### Abstract

A sketch-and-select Arnoldi process to generate a well-conditioned basis of a
Krylov space at low cost is proposed. At each iteration the procedure utilizes
randomized sketching to select a limited number of previously computed basis
vectors to project out of the current basis vector. The computational cost
grows linearly with the dimension of the Krylov space. The subset selection
problem for the projection step is approximately solved with a number of
heuristic algorithms and greedy methods used in statistical learning and
compressive sensing.

###### keywords:

Krylov method, Arnoldi process, randomized sketching

65F10, 65F50

##  1 Introduction

The Arnoldi process [2] is a key component of many Krylov subspace methods for
large-scale numerical linear algebra computations, including solving linear
systems of equations and eigenvalue problems with nonsymmetric matrices
$A\in\mathbb{R}^{N\times N}$; see, e.g., [30, 31]. The Arnoldi process is also
used for solving least squares problems, approximating matrix functions or
matrix equations, and in model order reduction, to name just a few other
applications. Given a starting vector ${\mathbf{b}}\in\mathbb{R}^{N}$ and an
integer $m\ll N$, the Arnoldi process iteratively constructs an orthonormal
basis $\\{{\mathbf{v}}_{1},{\mathbf{v}}_{2},\ldots,{\mathbf{v}}_{m}\\}$ of the
Krylov space
$\mathcal{K}_{m}(A,{\mathbf{b}}):=\operatorname{span}\\{{\mathbf{b}},A{\mathbf{b}},\ldots,A^{m-1}{\mathbf{b}}\\}$.
More precisely, given $j$ orthonormal basis vectors
${\mathbf{v}}_{1}:={\mathbf{b}}/\|{\mathbf{b}}\|,{\mathbf{v}}_{2},\ldots,{\mathbf{v}}_{j}$,
the next basis vector is obtained by orthogonalizing
${\mathbf{w}}_{j}:=A{\mathbf{v}}_{j}$ against all previous vectors,

(1)
$\widehat{\mathbf{w}}_{j}:={\mathbf{w}}_{j}-\sum_{i=1}^{j}h_{i,j}{\mathbf{v}}_{i},\quad
h_{i,j}:={\mathbf{v}}_{i}^{T}{\mathbf{w}}_{j},$

and then setting ${\mathbf{v}}_{j+1}:=\widehat{\mathbf{w}}_{j}/h_{j+1,j}$ with
$h_{j+1,j}=\|\widehat{\mathbf{w}}_{j}\|$. Collecting the basis vectors into
$V_{m}=[{\mathbf{v}}_{1},{\mathbf{v}}_{2},\ldots,{\mathbf{v}}_{m}]\in\mathbb{R}^{N\times
m}$ and the orthogonalization coefficients into
$H_{m}=[h_{i,j}]\in\mathbb{R}^{m\times m}$, the Arnoldi process generates an
Arnoldi decomposition

$AV_{m}=V_{m}H_{m}+h_{m+1,m}{\mathbf{v}}_{m+1}{\mathbf{e}}_{m}^{T},$

where ${\mathbf{e}}_{m}\in\mathbb{R}^{m}$ denotes the $m$-th canonical unit
vector. By construction, $H_{m}$ is an upper-Hessenberg matrix.

In terms of arithmetic cost, the Arnoldi process requires $m$ matrix-vector
products $A{\mathbf{w}}_{j}$, as well as $m(m+1)/2$ inner products and vector
operations (“axpy” in BLAS-1 naming), for a total of

$O(m\cdot\texttt{nnz}(A)+Nm^{2})$

arithmetic operations. For sufficiently sparse $A$, this cost will be
dominated by the $Nm^{2}$-term for the orthogonalization. There are at least
two possible ways to reduce this cost. The first one is to restart the Arnoldi
process after $m$ iterations, using ${\mathbf{v}}={\mathbf{v}}_{m+1}$ as the
starting vector for the next cycle. Such a restarting approach is particularly
natural in the context of solving linear systems of equations (as there exists
a linear error equation $A{\mathbf{e}}={\mathbf{r}}$ where ${\mathbf{r}}$ is
the residual), but it can also be used for eigenvalue problems [23, 35, 34] or
matrix function computations [1, 17]. Of course, the combined Krylov basis
computed after $\ell>1$ restarts is no longer orthonormal and this usually
leads to a delayed convergence in restarted Krylov methods.

The second, more recently proposed approach to reduce the arithmetic cost of
the Arnoldi process is to employ randomized sketching; see, e.g., [5, 6, 3,
24, 4]. The key tool of sketching is an embedding matrix
$S\in\mathbb{R}^{s\times N}$ with $m<s\ll N$ that distorts the Euclidean norm
$\|\,\cdot\,\|$ of vectors in a controlled manner [32, 38]. More precisely,
given a positive integer $m$ and some $\varepsilon\in[0,1)$, we assume that
$S$ is such that for all vectors ${\mathbf{v}}$ in the Krylov space
$\mathcal{K}_{m+1}(A,{\mathbf{b}})$,

(2)
$(1-\varepsilon)\|{\mathbf{v}}\|^{2}\leq\|S{\mathbf{v}}\|^{2}\leq(1+\varepsilon)\|{\mathbf{v}}\|^{2}.$

The matrix $S$ is called an _$\varepsilon$ -subspace embedding_ for
$\mathcal{K}_{m+1}(A,{\mathbf{b}})$. Condition Eq. 2 can equivalently be
stated with the Euclidean inner product [32, Cor. 4]: for all
${\mathbf{u}},{\mathbf{v}}\in\mathcal{K}_{m+1}(A,{\mathbf{b}})$,

${}\langle{\mathbf{u}},{\mathbf{v}}\rangle-\varepsilon\|{\mathbf{u}}\|\cdot\|{\mathbf{v}}\|\leq\langle
S{\mathbf{u}},S{\mathbf{v}}\rangle\leq\langle{\mathbf{u}},{\mathbf{v}}\rangle+\varepsilon\|{\mathbf{u}}\|\cdot\|{\mathbf{v}}\|.$

In practice, such a matrix $S$ is not explicitly available and we hence have
to draw it at random to achieve Eq. 2 with high probability. There are several
ways to construct a random matrix $S$ with this property, see for instance the
discussions in [24, Sec. 2.3] or [4, Sec. 2.1].

There are two main ways sketching can be employed within the Arnoldi process.
The first one, proposed and applied in [5, 6, 3], is to replace the inner
products computed in Eq. 1 by inner products on sketched vectors

$\widehat{\mathbf{w}}_{j}:={\mathbf{w}}_{j}-\sum_{i=1}^{j}\widetilde{h}_{i,j}{\mathbf{v}}_{i},\quad\widetilde{h}_{i,j}:=(S{\mathbf{v}}_{i})^{T}(S{\mathbf{w}}_{j}).$

Effectively, the process then computes an orthonormal sketched basis
$SV_{m+1}$. Using an efficient subspace embedding such as the subsampled
random cosine or Fourier transform [33, 39] requiring $O(N\log s)$ operations
when applied to a single vector, the overall complexity is now

$O(m\cdot\texttt{nnz}(A)+mN\log s+Nm^{2}).$

Even though the cost of computing all $(m+1)m/2$ inner products is reduced
from $N(m+1)m/2$ to $s(m+1)m/2$, there is still a quadratic dependency on $m$.
On the other hand, it follows from [4, Cor. 2.2] that

(3)
$\bigg{(}\frac{1-\varepsilon}{1+\varepsilon}\bigg{)}^{1/2}\operatorname{cond}(SV_{m+1})\leq\operatorname{cond}(V_{m+1})\leq\bigg{(}\frac{1+\varepsilon}{1-\varepsilon}\bigg{)}^{1/2}\operatorname{cond}(SV_{m+1}),$

so the computed Krylov basis $V_{m+1}$ will be close to orthonormal provided
that $\varepsilon$ is sufficiently small (i.e., $s$ is sufficiently large
which, in practice, means choosing $s$ between $2m$ and $4m$).

Alternatively, one may give up completely on the aim of computing a (near)
orthonormal Krylov basis and modify the target algorithm to deal with the non-
orthogonality gracefully. This has been proposed in [24] for sketched GMRES
and Rayleigh–Ritz extraction of eigenvalues, and in [19, 10] for matrix
function computations. One of the most straightforward approaches to generate
the Krylov basis with a reduced number of projection steps111Mathematically,
the term “orthogonalization” is no longer adequate when the basis is non-
orthogonal, so we use the term projection step to refer to a vector operation
in truncated Arnoldi. is the _truncated Arnoldi procedure._ Let a truncation
parameter $k$ be given, then in place of Eq. 1 we use the iteration

$\widehat{\mathbf{w}}_{j}:={\mathbf{w}}_{j}-\sum_{i=\max\\{1,j+1-k\\}}^{j}h_{i,j}{\mathbf{v}}_{i},\quad
h_{i,j}:={\mathbf{v}}_{i}^{T}{\mathbf{w}}_{j}.$

Alternatively, this same truncated Arnoldi procedure can be combined with
sketching by replacing the coefficients $h_{i,j}$ by their sketched
counterparts
$\widetilde{h}_{i,j}:=(S{\mathbf{v}}_{i})^{T}(S{\mathbf{w}}_{j})$. A key
benefit of truncation is that the $O(Nm^{2})$ cost of the orthogonalization is
reduced to $O(Nmk)$, i.e., it grows linearly with the Krylov basis dimension
$m$.

The truncated Arnoldi procedure is likely inspired by the Lanczos process for
a symmetric matrix $A$, which is mathematically equivalent to truncated
Arnoldi with $k=2$. The Faber–Manteuffel theorem [15] gives a complete
characterization of the matrices $A$ for which there is a short-term recursion
that generates an orthogonal set of Krylov basis vectors. The truncated
Arnoldi procedure first appeared in the context of eigenvalue problems [28,
Sec. 3.2] and linear systems [29, Sec. 3.3]. For computations involving matrix
functions, truncated Arnoldi has been used for the matrix exponential in [18,
20] and for more general matrix functions in [19, 10]. We also refer the
reader to [24, Sec. 4] and [27, Sec. 5].

Unfortunately, the Krylov bases generated by truncated Arnoldi (with and
without sketching) can become severely ill-conditioned even for moderate $m$.
In this paper we propose and test an alternative approach which we call the
_sketch-and-select Arnoldi process._ The key idea is simple: instead of
projecting each new Krylov basis vector against the $k$ previous basis
vectors, we use the sketched version of the Krylov basis to identify $k$
candidates for the projection. We show that this problem is related to a
sparse least squares approximation problem that has been studied in
statistical learning and compressive sensing. We then demonstrate with
performance profiles that the sketch-and-select Arnoldi process with a simple
select strategy to determine the candidate vectors often leads to much better
conditioned Krylov bases, outperforming truncated Arnoldi and many other
tested methods.

##  2 The sketch-and-select Arnoldi process

At iteration $j$ of the sketch-and-select Arnoldi process we compute
${\mathbf{w}}_{j}:=A{\mathbf{v}}_{j}$ as in the standard Arnoldi process, but
then aim to determine coefficients $h_{i,j}$ ($i=1,\ldots,j$) of which at most
$k$ are nonzero. We then compute

(4) $\widehat{\mathbf{w}}_{j}:={\mathbf{w}}_{j}-\sum_{h_{i,j}\neq
0}h_{i,j}{\mathbf{v}}_{i}$

and set ${\mathbf{v}}_{j+1}:=\widehat{\mathbf{w}}_{j}/h_{j+1,j}$ with a
suitable scaling factor $h_{j+1,j}$. To obtain the nonzero coefficients
$h_{i,j}$ we use the sketched Krylov basis $SV_{j}$ and the sketched vector
$S{\mathbf{w}}_{j}$. For any iteration $j>k$, we approximately solve the
following sparse least squares problem: select an index set
$I\subseteq\\{1,2,\ldots,j\\}$ with $|I|=k$ as

$\underset{I}{\arg\min}\underset{{\mathbf{h}}\in\mathbb{R}^{|I|}}{\min}\|S{\mathbf{w}}_{j}-SV_{j}(\,:\,,I){\mathbf{h}}\|.$

Here we have used MATLAB notation to denote column selection. Given the index
set $I$, the components of the best ${\mathbf{h}}$ are then used as the
projection coefficients $h_{i,j}$ in Eq. 4. Finally, the scaling factor
$h_{j+1,j}$ is chosen so that $\|S{\mathbf{v}}_{j+1}\|=1$.

The determination of an optimal index set $I$ is also known as best subset
selection problem, a classic topic of model selection in statistical learning
[22]. There are two main variants of this problem; (i) the problem considered
above where the sparsity level $k$ is prescribed, and (ii) for a given
tolerance $\epsilon>0$, the problem of selecting an index set $I$ with fewest
possible elements so that
$\min_{{\mathbf{h}}\in\mathbb{R}^{|I|}}\|S{\mathbf{w}}_{j}-SV_{j}(\,:\,,I){\mathbf{h}}\|\leq\epsilon$.
It is known that the determination of a global minimizer for such problems is
NP-hard; see [25]. Nevertheless, a vast amount of literature has been devoted
to developing efficient optimization algorithms for this task. A review of
these methods is beyond the scope of this paper, so we refer to the excellent
overview in [7].

A simple approach to selecting the index set $I$ is to retain the $k$
components of ${\mathbf{h}}:=(SV_{j})^{\dagger}(S{\mathbf{w}}_{j})$ which are
largest in modulus, ignoring the remaining $j-k$ components. We found this to
perform very well in the experiments reported in Section 3, and we present a
basic MATLAB implementation in Fig. 1. We refer to this variant as sketch +
select pinv (for “pseudoinverse”).

One could also try to justify a weighted version of sketch + select pinv by
noting that

$(SV_{j})^{T}[S{\mathbf{w}}_{j}-SV_{j}(\,:\,,I){\mathbf{h}}(I)-SV_{j}(\,:\,,\overline{I}){\mathbf{h}}(\overline{I})]=0,$

where $\overline{I}$ denotes the set of indices that have not been selected
and ${\mathbf{h}}=(SV_{j})^{\dagger}(S{\mathbf{w}}_{j})$ as before. Therefore

$(SV_{j})^{T}[\underbrace{S{\mathbf{w}}_{j}-SV_{j}(\,:\,,I){\mathbf{h}}(I)}_{=\,h_{j+1,j}S{\mathbf{v}}_{j+1}}]=(SV_{j})^{T}[SV_{j}(\,:\,,\overline{I}){\mathbf{h}}(\overline{I})].$

Our analysis in Section 2.1 shows that
$\operatorname{cond}([SV_{j},S{\mathbf{v}}_{j+1}])$ crucially depends on the
norm $\|(SV_{j})^{T}(S{\mathbf{v}}_{j+1})\|$, which should be kept as small as
possible. Upon applying norms we obtain

$|h_{j+1,j}|\cdot\|(SV_{j})^{T}(S{\mathbf{v}}_{j+1})\|\leq\sum_{i\in\overline{I}}|h_{i,j}|\cdot\|(SV_{j})^{T}(S{\mathbf{v}}_{i})\|,$

which suggests to assign to $I$ the indices of the $k$ coefficients $h_{i,j}$
for which $|h_{i,j}|\cdot\|(SV_{j})^{T}(S{\mathbf{v}}_{i})\|$ is largest
$(i=1,\ldots,j)$. Unfortunately, the non-orthogonality of $SV_{j}$ and the
complicated dependence of $h_{j+1,j}$ on the index set $I$ means that
minimizing the right-hand side of this inequality does not necessarily lead to
a significant reduction of $\|(SV_{j})^{T}(S{\mathbf{v}}_{j+1})\|$, which is
why we no longer consider this weighted variant of sketch + select pinv here.

Another approach is to select $I$ as above, but then to recompute the
corresponding coefficients as
${\mathbf{h}}=SV_{j}(\,:\,,I)^{\dagger}(S{\mathbf{w}}_{j})$. This ensures that
the projected $S\widehat{\mathbf{w}}_{j}$ is orthogonal to
$\operatorname{span}(SV_{j}(\,:\,,I))$. We refer to this variant as sketch +
select pinv2 (as two pseudoinverses are computed).

In the variant sketch + select corr (for “correlation”) we select $I$ as the
components of $(SV_{j})^{T}(S{\mathbf{w}}_{j})$ which are largest in modulus,
using the $k$ inner products as the projection coefficients $h_{i,j}$.222The
sketch + select corr variant has been hinted at in [24, version 1, Sec. 4.3]:
“Indeed,
$(\mathbf{S}\mathbf{q}_{i})^{*}(\mathbf{\mathbf{missing}}S\mathbf{A}\mathbf{q}_{j})\approx\mathbf{q}_{i}^{*}\mathbf{A}\mathbf{q}_{j}$
for all $i\leq j$, so we can choose to orthogonalize
$\mathbf{A}\mathbf{q}_{j}$ only against the basis vectors $\mathbf{q}_{i}$
where the inner product is nonnegligible.” Alternatively, we can recompute the
projection coefficients as
${\mathbf{h}}=SV_{j}(\,:\,,I)^{\dagger}(S{\mathbf{w}}_{j})$, referred to as
variant sketch + select corr pinv.

Finally, we also test three popular methods for sparse approximation, namely
orthogonal matching pursuit (OMP) [26, 37], subspace pursuit (SP) [11], and
the “Algorithm Greedy” from [25]. We have chosen these methods due to their
popularity but also because they naturally allow for a fixed sparsity level
$k$, as opposed to, e.g., LASSO [36].

function [V, H, SV, SAV] = ssa(A, b, m, k) %% Sketch-and-select Arnoldi
process sw = sketch(b); nsw = norm(sw); SV(:,1) = sw/nsw; V(:,1) = b/nsw; H =
[]; for j = 1:m w = A*V(:,j); sw = sketch(w); SAV(:,j) = sw; % the following
two lines perform the select operation coeffs = SV(:,1:j) \ sw; [coeffs,ind] =
maxk(abs(coeffs),k); w = w - V(:,ind)*coeffs; sw = sw - SV(:,ind)*coeffs; nsw
= norm(sw); SV(:,j+1) = sw/nsw; V(:,j+1) = w/nsw; H(ind,j) = coeffs; H(j+1,j)
= nsw; end

Figure 1: Basic (non-optimized) MATLAB implementation of the sketch-and-select
Arnoldi process. The function uses a sketch function that takes as input a
vector with $N$ components and returns a sketch with $s\ll N$ components. In
this variant, which we refer to as sketch + select pinv, the indices of Krylov
basis vectors to project out are obtained by keeping the $k$ coefficients of
the orthogonal projection which are largest in modulus. In an efficient
implementation, the least squares problems with the sketched basis $SV_{j}$
should be solved by retaining and updating at each iteration a QR
factorization of $SV_{j}$.

###  2.1 Growth of the basis condition number

By Eq. 3 it is sufficient to only control the condition number growth of the
sketched Krylov basis $[SV_{j},S{\mathbf{v}}_{j+1}]$. For notational
convenience, we will write this as $[V,{\mathbf{v}}]$. Our aim is to bound the
growth of $\operatorname{cond}([V,{\mathbf{v}}])$ in terms of the current
$\operatorname{cond}(V)$. Note that the Gram matrix
$[V,{\mathbf{v}}]^{T}[V,{\mathbf{v}}]$ has unit diagonal in the sketch-and-
select Arnoldi process as we always normalize each sketched Krylov basis
vector, and so

$\operatorname{cond}(V)^{2}=\operatorname{cond}\left(\begin{bmatrix}V^{T}V&0\\\
0^{T}&1\end{bmatrix}\right).$

Also,

$\operatorname{cond}([V,{\mathbf{v}}])^{2}=\operatorname{cond}\left(\begin{bmatrix}V^{T}V&V^{T}{\mathbf{v}}\\\
{\mathbf{v}}^{T}V&1\end{bmatrix}\right).$

We can apply standard relative perturbation bounds known for symmetric
positive definite matrices. To this end, write

$\begin{bmatrix}V^{T}V&V^{T}{\mathbf{v}}\\\
{\mathbf{v}}^{T}V&1\end{bmatrix}=\begin{bmatrix}V^{T}V&0\\\
0^{T}&1\end{bmatrix}+\begin{bmatrix}O&V^{T}{\mathbf{v}}\\\
{\mathbf{v}}^{T}V&0\end{bmatrix}=:G+\Delta G.$

Note that $\|\Delta G\|=\|V^{T}{\mathbf{v}}\|$. Then (see, e.g., [13] or [21,
Thm. 2.3])

$\displaystyle\operatorname{cond}([V,{\mathbf{v}}])^{2}$ $\displaystyle=$
$\displaystyle\operatorname{cond}(G+\Delta G)=\frac{\lambda_{\max}(G+\Delta
G)}{\lambda_{\min}(G+\Delta G)}$ $\displaystyle\leq$
$\displaystyle\frac{(1+\eta)\lambda_{\max}(G)}{(1-\eta)\lambda_{\min}(G)}=\frac{1+\eta}{1-\eta}\operatorname{cond}(V)^{2},$

where

$\eta=\|G^{-1/2}(\Delta
G)G^{-1/2}\|=\|(V^{T}V)^{-1/2}V^{T}{\mathbf{v}}\|\leq\sigma_{\min}(V)^{-1}\|V^{T}{\mathbf{v}}\|.$

Clearly this bound is only useful as long as $1-\eta>0$, which is guaranteed
if $\|V^{T}{\mathbf{v}}\|<\sigma_{\min}(V)$, and generally it cannot be
expected to be sharp. However, it shows that it is a good idea to try to keep
$\|V^{T}{\mathbf{v}}\|$ as small as possible. Going back to the notation used
for the sketch-and-select Arnoldi process, this means that we should aim to
keep $\|(SV_{j})^{T}(S{\mathbf{v}}_{j+1})\|$ small.

An alternative approach to quantify the condition number growth is by bounding
the smallest and largest eigenvalue of the Gram matrix, taking into account
the special structure of that matrix. The following theorem provides such a
result, giving a more explicit bound on
$\operatorname{cond}([V,{\mathbf{v}}])$.

######  Theorem 2.1.

Let $V$ be a matrix with $m$ linearly independent columns of unit norm. Denote
by $\sigma_{\min}$ and $\sigma_{\max}$ the smallest and largest singular value
of $V$, respectively. Further, let ${\mathbf{v}}$ be a unit norm vector such
that $\|V^{T}{\mathbf{v}}\|<\sigma_{\min}$. Then

(5)
$\operatorname{cond}([V,{\mathbf{v}}])^{2}\leq\frac{1+{\sigma_{\max}^{2}}+\sqrt{(\sigma_{\max}^{2}-1)^{2}+4\|V^{T}{\mathbf{v}}\|^{2}}}{1+\sigma_{\min}^{2}-\sqrt{(\sigma_{\min}^{2}-1)^{2}+4\|V^{T}{\mathbf{v}}\|^{2}}}.$

######  Proof 2.2.

We begin by remarking that, as the columns of $[V,{\mathbf{v}}]$ are
normalized, we have $\sigma_{\max}([V,{\mathbf{v}}])\leq\sqrt{m+1}$ and hence
any ill-conditioning of $[V,{\mathbf{v}}]$ is mainly attributable to a small
value of

$\sigma_{\min}([V,{\mathbf{v}}])=\lambda_{\min}([V,{\mathbf{v}}]^{T}[V,{\mathbf{v}}])^{1/2}.$

Define the Rayleigh quotient

$R({\mathbf{x}},{\mathbf{v}})=\left[{\mathbf{x}}^{T},\sqrt{1-\|{\mathbf{x}}\|^{2}}\right]\begin{bmatrix}V^{T}V&V^{T}{\mathbf{v}}\\\
{\mathbf{v}}^{T}V&1\end{bmatrix}\begin{bmatrix}{\mathbf{x}}\\\
\sqrt{1-\|{\mathbf{x}}\|^{2}}\end{bmatrix},\quad\|{\mathbf{x}}\|\leq 1.$

Let us denote $\beta=\|{\mathbf{x}}\|$ and fix
$\|V^{T}{\mathbf{v}}\|:=\alpha$. Then

(6)
$R({\mathbf{x}},{\mathbf{v}})={\mathbf{x}}^{T}(V^{T}V){\mathbf{x}}+2{\mathbf{x}}^{T}(V^{T}{\mathbf{v}})\sqrt{1-\beta^{2}}+1-\beta^{2}.$

The first term in this Rayleigh quotient is minimized by choosing
${\mathbf{x}}$ as an eigenvector of $V^{T}V$ corresponding to
$\lambda_{\min}=\lambda_{\min}(V^{T}V)$. For any choice of ${\mathbf{x}}$,
$\|{\mathbf{x}}\|=\beta$, the second term is minimal when ${\mathbf{v}}$ is
such that $V^{T}{\mathbf{v}}=-\alpha{\mathbf{x}}/\beta$. Hence, we can
minimize the overall Rayleigh quotient directly, leading to

$R_{\min}(\beta)=\beta^{2}\lambda_{\min}-2\alpha\beta\sqrt{1-\beta^{2}}+1-\beta^{2}.$

To find the optimal $\beta\in[0,1]$, we set $\gamma:=1-\beta^{2}$, upon which
$R_{\min}(\beta)=(1-\gamma)\lambda_{\min}-2\alpha\sqrt{\gamma-\gamma^{2}}+\gamma$
which is easy to differentiate for $\gamma$. The optimal value $\beta_{*}$
that minimizes $R_{\min}$ is

$\beta_{*}=\sqrt{1-\gamma_{*}},\quad\text{where}\quad\gamma_{*}=\frac{C_{*}^{2}-C_{*}\sqrt{C_{*}^{2}+4}+4}{2(C_{*}^{2}+4)},\quad
C_{*}=\frac{1-\lambda_{\min}}{\alpha}.$

We can now derive a rather simple expression for $R_{\text{min}}(\beta_{*})$
in terms of $C_{*}$. We have

$\gamma_{*}=\frac{C_{*}^{2}-C_{*}\sqrt{C_{*}^{2}+4}+4}{2(C_{*}^{2}+4)}=\frac{1}{2}-\frac{1}{2}\frac{C_{*}}{\sqrt{C_{*}^{2}+4}},$

and from this expression it easily follows that

$\sqrt{\gamma_{*}-\gamma_{*}^{2}}=\frac{1}{\sqrt{C_{*}^{2}+4}}.$

By plugging this expression and $\gamma_{*}$ into the expression of
$R_{\text{min}}(\beta_{*})$, we have

$\displaystyle R_{\text{min}}(\beta_{*})$
$\displaystyle=\lambda_{\text{min}}+(1-\lambda_{\text{min}})\gamma_{*}-2\alpha\sqrt{\gamma_{*}-\gamma_{*}^{2}}$
$\displaystyle=\lambda_{\text{min}}+\alpha
C_{*}\Big{(}\frac{1}{2}-\frac{1}{2}\frac{C_{*}}{\sqrt{C_{*}^{2}+4}}\Big{)}-\frac{2\alpha}{\sqrt{C_{*}^{2}+4}}$
$\displaystyle=\lambda_{\text{min}}+\frac{\alpha}{2\sqrt{C_{*}^{2}+4}}\Big{(}C_{*}\sqrt{C_{*}^{2}+4}-C_{*}^{2}-4\Big{)}$
$\displaystyle=\lambda_{\text{min}}+\frac{\alpha}{2}\Big{(}C_{*}-\sqrt{C_{*}^{2}+4}\Big{)}$
$\displaystyle=\frac{1}{2}+\frac{1}{2}\lambda_{\text{min}}-\frac{\alpha}{2}\sqrt{C_{*}^{2}+4}$
$\displaystyle=\frac{1+\lambda_{\text{min}}-\sqrt{(1-\lambda_{\text{min}})^{2}+4\alpha^{2}}}{2}.$

For this quantity to be positive, we require $\alpha^{2}<\lambda_{\min}$ or
equivalently $\|V^{T}v\|<\sigma_{\min}(V)$.

Similarly, the first term in the Rayleigh quotient Eq. 6 is _maximized_ by
choosing ${\mathbf{x}}$ as an eigenvector of $V^{T}V$ corresponding to
$\lambda_{\max}=\lambda_{\max}(V^{T}V)$. For any choice of ${\mathbf{x}}$,
$\|{\mathbf{x}}\|=\beta$, the second term in Eq. 6 is maximal when
${\mathbf{v}}$ is such that $V^{T}{\mathbf{v}}=\alpha{\mathbf{x}}/\beta$.
Hence, we can maximize the overall Rayleigh quotient directly, leading to

$R_{\max}(\beta)=\beta^{2}\lambda_{\max}+2\alpha\beta\sqrt{1-\beta^{2}}+1-\beta^{2}.$

To find the optimal $\beta\in[0,1]$, we set $\gamma:=1-\beta^{2}$, upon which
$R_{\max}(\beta)=(1-\gamma)\lambda_{\max}+2\alpha\sqrt{\gamma-\gamma^{2}}+\gamma$.
The optimal value $\beta^{*}$ that maximizes $R_{\max}$ is

$\beta^{*}=\sqrt{1-\gamma^{*}},\quad\text{where}\quad\gamma^{*}=\frac{C^{2}-C\sqrt{C^{2}+4}+4}{2(C^{2}+4)},\quad
C=\frac{\lambda_{\max}-1}{\alpha}.$

(Note that the only difference compared to the above is in $C$ versus
$C_{*}$.) Evaluating $R_{\max}(\beta^{*})$ yields

$R_{\text{max}}(\beta^{*})=\frac{1+\lambda_{\text{max}}+\sqrt{(\lambda_{\text{max}}-1)^{2}+4\alpha^{2}}}{2}.$

Combining the expressions for the (worst-case) Rayleigh quotients, we obtain

$\operatorname{cond}([V,{\mathbf{v}}])^{2}\leq\frac{R_{\max}(\beta^{*})}{R_{\min}(\beta_{*})}=\frac{1+\lambda_{\text{max}}+\sqrt{(\lambda_{\text{max}}-1)^{2}+4\alpha^{2}}}{1+\lambda_{\text{min}}-\sqrt{(\lambda_{\text{min}}-1)^{2}+4\alpha^{2}}}.$

The result follows since $\lambda_{\max}=\sigma_{\max}(V)^{2}$,
$\lambda_{\min}=\sigma_{\min}(V)^{2}$, and $\alpha=\|V^{T}{\mathbf{v}}\|$.

Observe that it follows from the proof of Theorem 2.1 that if the vector
${\mathbf{v}}$ minimizes the Rayleigh quotient in Eq. 6 (i.e., it satisfies
$V^{T}{\mathbf{v}}=-\alpha{\mathbf{x}}/\beta_{*}$ where ${\mathbf{x}}$ is the
eigenvector of $V^{T}V$ corresponding to $\lambda_{\text{min}}$), we have

(7)
$\sigma_{\min}([V,{\mathbf{v}}])^{2}=\frac{1+\sigma_{\text{min}}^{2}-\sqrt{(1-\sigma_{\text{min}}^{2})^{2}+4\|V^{T}{\mathbf{v}}\|^{2}}}{2}.$

Since for any vector ${\mathbf{v}}$ of unit norm we have
$\sigma_{\max}([V,{\mathbf{v}}])\geq\sigma_{\max}\geq 1$, this implies that
there is a choice of ${\mathbf{v}}$ for which

(8)
$\operatorname{cond}([V,{\mathbf{v}}])^{2}\geq\frac{2}{1+\sigma_{\min}^{2}-\sqrt{(\sigma_{\min}^{2}-1)^{2}+4\|V^{T}{\mathbf{v}}\|^{2}}}.$

Recalling that $\sigma_{\text{max}}([V,{\mathbf{v}}])\leq\sqrt{m+1}$, we find
that the right-hand side of Eq. 8 is smaller than the right-hand side of Eq. 5
at most by a factor $m+1$. In principle, it is possible to select a vector
${\mathbf{v}}$ that realizes Eq. 8 at every iteration, even though this is
unlikely to happen in practice.

Going back to the notation used for the sketch-and-select Arnoldi process, let
us now consider the behavior of $\sigma_{\text{min}}(SV_{m})$ as $m$
increases, assuming that at each iteration $S{\mathbf{v}}_{m+1}$ is selected
as a vector ${\mathbf{v}}$ that satisfies Eq. 7. For convenience, define
$x_{m}:=\sigma_{\text{min}}(SV_{m})$ and
$\alpha_{m}:=\|S{\mathbf{v}}_{m+1}\|<x_{m}$. Because of Eq. 7, these
quantities satisfy the recurrence relation

$x_{m+1}^{2}=\frac{1}{2}\left(1+x_{m}^{2}-\sqrt{(1-x_{m}^{2})^{2}+4\alpha_{m}^{2}}\right),\qquad
m\geq 1,$

with $x_{1}=\sigma_{\text{min}}(S{\mathbf{v}}_{1})=1$. Using the fact that
$\sqrt{1+z}\geq 1+\frac{1}{2}z-\frac{1}{8}z^{2}$ for all $z\geq 0$, we can
show that

$x_{m+1}^{2}\leq
x_{m}^{2}-\frac{\alpha_{m}^{2}}{1-x_{m}^{2}}\left(1-\frac{\alpha_{m}^{2}}{(1-x_{m}^{2})^{2}}\right).$

If for instance we take $x_{m}\leq 1/\sqrt{2}$ and
$\frac{1}{2}x_{m}\leq\alpha_{m}<x_{m}$, we have

$1-\frac{\alpha_{m}^{2}}{(1-x_{m}^{2})^{2}}\leq\frac{1}{2},$

and so we obtain

$x_{m+1}^{2}\leq
x_{m}^{2}-\frac{1}{2}\alpha_{m}^{2}\leq\frac{7}{8}x_{m}^{2},\qquad m\geq 1.$

This implies that

(9) $\sigma_{\text{min}}(SV_{m})\leq\left(\frac{7}{8}\right)^{m/2},\qquad
m\geq 1,$

showing that a geometric convergence to $0$ of the smallest singular value of
$SV_{m}$ is possible. As a consequence, the condition number of $SV_{m}$ may
diverge geometrically (and hence also $\operatorname{cond}(V_{m})$, because of
Eq. 3), even if we impose the condition
$\lVert(SV_{m})^{T}S{\mathbf{v}}_{m+1}\rVert=\frac{1}{2}\sigma_{\text{min}}(SV_{m})$,
which is quite stringent.

##  3 Numerical experiments

We now test variants of the proposed sketch-and-select Arnoldi process on a
range of matrices from the SuiteSparse Matrix Collection (formerly the
University of Florida Sparse Matrix Collection [12]), a widely used set of
sparse matrix benchmarks collected from a wide range of applications. We
include 80 matrices $A$ in our test, which correspond to all square
numerically nonsymmetric matrices in the collection (as of June 2023) with
sizes between $N=10^{4}$ and $10^{6}$. The starting vector ${\mathbf{b}}$ is
chosen at random with unit normally distributed entries and kept constant for
all tests with the same matrix dimension. The MATLAB scripts to reproduce the
experiments in this section are available at
https://github.com/simunec/sketch-select-arnoldi.

In the tests below we compare the seven variants of the sketch-and-select
Arnoldi process introduced in Section 2, each using a different method for the
sparse least squares problem, as well as the truncated Arnoldi process with
and without sketching.

###  3.1 Illustration with a single matrix

For our first experiment we plot in Fig. 2 the condition number
$\operatorname{cond}(V_{m})$ as a function of $m$ for the SuiteSparse problem
Norris/torso3, a matrix of size $N=259,156$ corresponding to a finite
difference electro-physiological 3D model of a torso. We use a truncation
parameter of $k=2$ and $k=5$ and perform $m=100$ Arnoldi iterations. The
sketching operator is the subsampled random Hadamard transform (SRHT) [3, Sec.
1.3] with an embedding dimension of $s=200$.

We see from Fig. 2 that most sketch-and-select Arnoldi variants exhibit a much
smaller condition number growth than truncated Arnoldi. The main exception is
sketch + select corr which performs very badly in both cases, and sketch +
select corr pinv which is only acceptable for $k=5$. Surprisingly, the variant
sketch + select pinv2 breaks down after about $m=80$ Arnoldi iterations in the
case $k=2$, but it works very well for $k=5$. The four most reliable variants
are sketch + select pinv, sketch + select OMP, sketch + select SP, and sketch
+ select greedy, all leading to a rather slow growth of the condition number.

In terms of computational cost, all sketch-and-select variants only perform
operations on small sketched matrices and vectors to determine the projection
coefficients. Hence these computations are comparably cheap, but sketch +
select pinv is the cheapest method as it only requires the solution of a
single $s\times j$ least squares problem in the $j$-th Arnoldi iteration. OMP
is a greedy method that, for each Arnoldi iteration $j$, requires $k$
iterations, with the $i$-th inner iteration involving a matrix-vector product
with the $s\times j$ matrix $SV_{j}$ and the solution of an $s\times i$ least
squares problem with the currently selected columns of $SV_{j}$
($i=1,\ldots,j$). SP is also an iterative method but we have fixed the number
of iterations to 1. As a result, the operations performed by SP in the $j$-th
Arnoldi iteration are two matrix-vector products with the $s\times j$ matrix
$SV_{j}$, a matrix-vector product with an $s\times k$ matrix, and the solution
of two $s\times k$ least squares problems, as well as a least squares problem
with a matrix of size at most $s\times 2k$ formed with selected columns of
$SV_{j}$. The “Algorithm Greedy” from [25] is an iterative method that
requires $k$ iterations for each Arnoldi iteration $j$, and each inner
iteration involves the computation of $2j$ inner products between sketched
vectors, for a total cost of $2jk$ sketched inner products and the solution of
one $s\times k$ least squares problem in the $j$-th Arnoldi iteration.

Figure 2: Basis condition number growth of nine different methods for the
matrix Norris/torso3, using a truncation parameter of $k=2$ (top) and $k=5$
(bottom).

###  3.2 Performance profiles

Our next tests involve all 80 matrices of the SuiteSparse Matrix Collection
and we use performance profiles [14] to visualize the results. In a
performance profile, each algorithm is represented by a non-decreasing curve
in a $\theta$–$y$ graph. The $\theta$-axis represents a tolerance $\theta\geq
1$ and the $y$-axis corresponds to a fraction $y\in[0,1]$. If a curve passes
through a point $(\theta,y)$ this means that the corresponding algorithm
performed within a factor $\theta$ of the best observed performance on $100y$%
of the test problems. For $\theta=1$ one can read off on what fraction of all
test problems each algorithm was the best performer, while as
$\theta\to\infty$ all curves approach the value $y=1$, unless an algorithm has
failed on a fraction of the test problems.

For each test problem, we run each algorithm until the condition number of the
constructed basis becomes larger than $10^{12}$, up to a maximum basis
dimension. The performance ratio is computed as the inverse of the basis
dimension that is reached, so that, e.g., $\theta=2$ would correspond to an
algorithm that generates a Krylov basis with half the size of the basis
generated by the best algorithm.

The top panel in Fig. 3 shows the performance profiles for a target basis
condition number of $10^{12}$ and a maximum basis dimension of $m=100$, with a
truncation parameter of $k=2$. The embedding dimension is $s=200$. The four
most reliable variants are sketch + select pinv, sketch + select OMP, sketch +
select SP, and sketch + select greedy, and they are almost indistinguishable
in performance. Given that sketch + select pinv is the most straightforward to
implement and the most computationally efficient, it emerges as the method of
choice from these experiments. The bottom panel in Fig. 3 displays the
dimensions of the Krylov bases constructed for each test matrix by the four
best performing algorithms and by truncated Arnoldi. The matrices are sorted
so that the dimensions of the bases generated by truncated Arnoldi are in non-
decreasing order.

Figure 3: Basis sizes reached by the different methods, using a truncation
parameter of $k=2$, with maximum basis size $m=100$ and condition number
bounded by $10^{12}$. Performance profiles (top) and basis sizes for the best
performing methods (bottom).

A similar picture emerges in Figs. 4 and 5, where we have increased the
truncation parameter to $k=5$ and $k=10$, respectively; the maximum basis
dimension is increased to $m=150$ and $m=200$, respectively, and the embedding
dimension is chosen as $s=2m$. With these parameters, the difference in
performance between truncated Arnoldi and the best sketch + select variants
becomes more significant, and the sketch + select pinv variant can be seen to
perform slightly better than sketch + select OMP, sketch + select SP and
sketch + select greedy.

Figure 4: Basis sizes reached by the different methods, using a truncation
parameter of $k=5$, with maximum basis size $m=150$ and condition number
bounded by $10^{12}$. Performance profiles (top) and basis sizes for the best
performing methods (bottom).

Figure 5: Basis sizes reached by the different methods, using a truncation
parameter of $k=10$, with maximum basis size $m=200$ and condition number
bounded by $10^{12}$. Performance profiles (top) and basis sizes for the best
performing methods (bottom).

###  3.3 The effect of the starting vector

In the previous experiments, we have used a starting vector $\boldsymbol{b}$
with random unit normally distributed entries. To investigate the influence of
the starting vector on the performance of the algorithms, we repeat the
experiment of Fig. 4 using as vector ${\mathbf{b}}$ the first canonical unit
vector ${\mathbf{e}}_{1}$, instead of a random vector. The results are
reported in Fig. 6. Surprisingly, the sketch-and-select Arnoldi variants
perform significantly worse with this starting vector while, on the other
hand, the overall performance of truncated Arnoldi is almost the same as
before.

In Fig. 7 we repeat the same experiment by slightly perturbing the starting
vector, i.e., we take ${\mathbf{b}}={\mathbf{e}}_{1}+10^{-15}{\mathbf{e}}$,
where ${\mathbf{e}}$ denotes the vector of all ones. With this change, the
performance of the sketch-and-select Arnoldi variants improves significantly
relative to truncated Arnoldi, though not to the same level of what was
observed in the experiment in Fig. 4 with a random starting vector. A very
similar improvement is obtained with a small random perturbation of the vector
${\mathbf{b}}$.

While it currently appears to be impossible to make any general statements
about the dependence of relative performance of sketch-and-select Arnoldi on
the starting vector ${\mathbf{b}}$, we observed that truncated Arnoldi can
produce artificially well conditioned bases for certain starting vectors. For
example, it may happen that sparse basis vectors constructed by truncated
Arnoldi have disjoint supports, and so they are all orthogonal to each other.
(One example is the matrix Norris/lung2: when ${\mathbf{b}}={\mathbf{e}}_{1}$,
the first 452 Krylov basis vectors produced by truncated Arnoldi with
truncation parameter $k\geq 2$ are given by $\pm{\mathbf{e}}_{2j-1}$.) The
sketching-based methods do not “see” the sparsity of the basis vectors and
hence cannot produce this exact orthogonality, losing performance relative to
truncated Arnoldi. Adding a small (random) perturbation to the starting vector
${\mathbf{b}}$ removes the sparsity and hence reduces the appearance of such
artificial cases.

Figure 6: Basis sizes reached by the different methods, using a truncation
parameter of $k=5$, with maximum basis size $m=150$ and condition number
bounded by $10^{12}$. The starting vector is ${\mathbf{b}}={\mathbf{e}}_{1}$.
Performance profiles (top) and basis sizes for the best performing methods
(bottom).

Figure 7: Basis sizes reached by the different methods, using a truncation
parameter of $k=10$, with maximum basis size $m=150$ and condition number
bounded by $10^{12}$. The starting vector is
${\mathbf{b}}={\mathbf{e}}_{1}+10^{-15}{\mathbf{e}}$. Performance profiles
(top) and basis sizes for the best performing methods (bottom).

##  4 Further remarks on the subset selection problem

Selecting $k$ columns that give the smallest condition number of
$[V,{\mathbf{v}}]$ in the sketch-and-select Arnoldi process is a combinatorial
problem, and even selecting a near-best index set $I$ is nontrivial. We now
give examples demonstrating that neither the largest coefficients in
$V^{\dagger}{\mathbf{w}}$ nor those in $V^{T}{\mathbf{w}}$ do necessarily
indicate the best vectors to select for a minimal condition number growth.
Note that most greedy algorithms, including the “Algorithm Greedy” in [25],
OMP [26], and SP [11] use the entries of either $V^{\dagger}{\mathbf{w}}$ or
$V^{T}{\mathbf{w}}$ to select vectors, and so can be misled on examples like
the ones below. Consider

$V=\frac{1}{\sqrt{5}}\begin{bmatrix}1&0&0\\\ 2&2&0\\\ 0&1&1\\\
0&0&2\end{bmatrix},\quad\boldsymbol{w}=\begin{bmatrix}8\\\ 8\\\ 9\\\
7\end{bmatrix}.$

Note that all columns of $V$ have unit norm. Say $k=1$, then which of the
three columns of $V$ should we project out of ${\mathbf{w}}$ to get the best
possible conditioned basis of four vectors? Let us compute

$V^{\dagger}\boldsymbol{w}=\begin{bmatrix}9.39\\\ 1.68\\\
9.95\end{bmatrix},\quad V^{T}{\mathbf{w}}=\begin{bmatrix}10.7\\\ 11.2\\\
10.3\end{bmatrix}.$

According to these coefficients, we might be tempted to project ${\mathbf{w}}$
either against $\boldsymbol{v}_{3}$ or $\boldsymbol{v}_{2}$, respectively.
However, among the vectors

$\widehat{\boldsymbol{v}}_{i}:=(I-\boldsymbol{v}_{i}\boldsymbol{v}_{i}^{\dagger}){\mathbf{w}},\quad
i=1,2,3,$

it is the choice $i=1$ that strictly minimizes both
$\operatorname{cond}([V,\widehat{\boldsymbol{v}}_{i}])$ and
$\operatorname{cond}([V,\widehat{\boldsymbol{v}}_{i}/{\|\widehat{\boldsymbol{v}}_{i}\|}])$.

Instead of condition number growth, perhaps a better measure to look at is the
growth of loss of orthogonality, e.g., by any of the metrics

$\left\|I-[V,\boldsymbol{v}]^{T}[V,\boldsymbol{v}]\right\|,\quad\left\|I-\left[V,\frac{\boldsymbol{v}}{\|\boldsymbol{v}\|}\right]^{T}\left[V,\frac{\boldsymbol{v}}{\|\boldsymbol{v}\|}\right]\right\|,$
$\left\|[V,\boldsymbol{v}]^{T}[V,\boldsymbol{v}]\right\|,\quad\left\|\left[V,\frac{\boldsymbol{v}}{\|\boldsymbol{v}\|}\right]^{T}\left[V,\frac{\boldsymbol{v}}{\|\boldsymbol{v}\|}\right]\right\|$

in the Euclidean or Frobenius norm. Can we get some guarantees for that? It
turns out that this is also not the case. The above matrix $V$ and the vector
${\mathbf{w}}=[9,9,10,10]^{T}$ give examples where, in all cases, the third
component of $V^{\dagger}\boldsymbol{w}$ and $V^{T}\boldsymbol{w}$ is the
largest in modulus, but the smallest growth in loss of orthogonality with
$k=1$ projection steps is obtained by projecting ${\mathbf{w}}$ against
${\mathbf{v}}_{2}$.

We conclude that there must be some condition on $V$ (and possibly
${\mathbf{w}}$) to guarantee that the “correct” (optimal) support of $k$
coefficients is selected. This is a well-known fact in compressive sensing
[16], where conditions like the restricted isometry property (RIP, [9]) are
needed to guarantee exact or approximate recovery in sparse least squares
approximation (see, e.g., [8] or [37, Thm. C]). On the other hand, the sensing
problem $A{\mathbf{x}}={\mathbf{b}}$ studied in this field is usually
underdetermined and one often has the freedom to choose the columns of $A$
(the dictionary), whereas in our case we would like to select $k$ columns from
a basis which is otherwise unstructured. Moreover, in compressive sensing the
vector ${\mathbf{x}}$ that one wants to recover is usually sparse, or at least
well approximated by a sparse vector, while in our case the coefficient
vectors are generally dense. Nevertheless, we believe that results developed
in compressive sensing may be used to gain more insight into the sketch-and-
select Arnoldi process.

##  5 Conclusions and future work

We have introduced a sketch-and-select Arnoldi process and demonstrated its
potential to generate Krylov bases that are significantly better conditioned
than those computed with the truncated Arnoldi process, at a computational
cost that grows only linearly with the dimension of the Krylov space. We have
identified that the problem of generating a well-conditioned Krylov basis in
that way is related to the best subset selection problem in statistical
learning and the sparse approximation problem encountered in compressive
sensing. While in principle any of the many methods that have been proposed
for these problems can be employed, we have been surprised to find that the
most basic variant of the sketch-and-select Arnoldi process shown in Fig. 1 is
among the best. In this approach we simply retain the $k$ largest modulus
coefficients of the least squares solution, setting the remaining coefficients
to zero.

Our implementation of the sketch-and-select Arnoldi process in Fig. 1 is not
optimized for performance, with several straightforward improvements possible
including the use of QR updating strategies for the solution of the least
squares problem or performing some of the operations in reduced precision.
Also, for a practical implementation to be used in production, the sketch-and-
select Arnoldi process could be modified to adapt the parameter $k$
dynamically based on the measured growth of the condition number. Bounds like
the ones derived in Section 2.1 might be useful to control the condition
number growth efficiently. Further possible extensions include a sketch-and-
select block Arnoldi process or restarting strategies.

## Acknowledgments

I. S. wishes to thank S. G. for the kind hospitality received during a
research visit to the University of Manchester, where this work was initiated.

## References

  * [1] Martin Afanasjew, Michael Eiermann, Oliver G. Ernst, and Stefan Güttel.  Implementation of a restarted Krylov subspace method for the evaluation of matrix functions.  Linear Algebra Appl., 429:2293–2314, 2008.
  * [2] Walter E. Arnoldi.  The principle of minimized iteration in the solution of the matrix eigenvalue problem.  Q. Appl. Math., 9:17–29, 1951.
  * [3] Oleg Balabanov and Laura Grigori.  Randomized block Gram–Schmidt process for solution of linear systems and eigenvalue problems.  arXiv preprint arXiv:2111.14641, 2021.
  * [4] Oleg Balabanov and Laura Grigori.  Randomized Gram–Schmidt process with application to GMRES.  SIAM J. Sci. Comput., 44(3):A1450–A1474, 2022.
  * [5] Oleg Balabanov and Anthony Nouy.  Randomized linear algebra for model reduction. Part I: Galerkin methods and error estimation.  Adv. Comput. Math., 45:2969–3019, 2019.
  * [6] Oleg Balabanov and Anthony Nouy.  Randomized linear algebra for model reduction–part II: minimal residual methods and dictionary-based approximation.  Adv. Comput. Math., 47:1–54, 2021.
  * [7] Dimitris Bertsimas, Angela King, and Rahul Mazumder.  Best subset selection via a modern optimization lens.  Ann. Statist., 44(2):813–852, 2016.
  * [8] Emmanuel J. Candes, Justin K. Romberg, and Terence Tao.  Stable signal recovery from incomplete and inaccurate measurements.  Comm. Pure Appl. Math., 59(8):1207–1223, 2006.
  * [9] Emmanuel J. Candes and Terence Tao.  Decoding by linear programming.  IEEE Trans. Inform. Theory, 51(12):4203–4215, 2005.
  * [10] Alice Cortinovis, Daniel Kressner, and Yuji Nakatsukasa.  Speeding up Krylov subspace methods for computing $f({A})b$ via randomization.  arXiv preprint arXiv:2212.12758, 2022.
  * [11] Wei Dai and Olgica Milenkovic.  Subspace pursuit for compressive sensing signal reconstruction.  IEEE Trans. Inform. Theory, 55(5):2230–2249, 2009.
  * [12] Timothy A. Davis and Yifan Hu.  The University of Florida sparse matrix collection.  ACM Trans. Math. Software, 38(1):1–25, 2011.
  * [13] James Demmel and Krešimir Veselić.  Jacobi’s method is more accurate than QR.  SIAM J. Matrix Anal. Appl., 13(4):1204–1245, 1992.
  * [14] Elizabeth D. Dolan and Jorge J. Moré.  Benchmarking optimization software with performance profiles.  Math. Program., 91:201–213, 2002.
  * [15] Vance Faber and Thomas Manteuffel.  Necessary and sufficient conditions for the existence of a conjugate gradient method.  SIAM J. Numer. Anal., 21(2):352–362, 1984.
  * [16] Simon Foucart and Holger Rauhut.  An Invitation to Compressive Sensing.  Springer, 2013.
  * [17] Andreas Frommer, Stefan Güttel, and Marcel Schweitzer.  Efficient and stable Arnoldi restarts for matrix functions based on quadrature.  SIAM J. Matrix Anal. Appl., 35:661–683, 2014.
  * [18] Stéphane Gaudreault, Greg Rainwater, and Mayya Tokman.  KIOPS: a fast adaptive Krylov subspace solver for exponential integrators.  J. Comput. Phys., 372:236–255, 2018.
  * [19] Stefan Güttel and Marcel Schweitzer.  Randomized sketching for Krylov approximations of large-scale matrix functions.  To appear in SIAM J. Matrix Anal. Appl., 2023.
  * [20] Antti Koskela.  Approximating the matrix exponential of an advection-diffusion operator using the incomplete orthogonalization method.  In Numerical Mathematics and Advanced Applications-ENUMATH 2013, pages 345–353. Springer, 2015.
  * [21] Roy Mathias.  Spectral perturbation bounds for positive definite matrices.  SIAM J. Matrix Anal. Appl., 18(4):959–980, 1997.
  * [22] Alan Miller.  Subset Selection in Regression.  CRC Press, 2002.
  * [23] Ronald B. Morgan.  On restarting the Arnoldi method for large nonsymmetric eigenvalue problems.  Math. Comp., 65(215):1213–1230, 1996.
  * [24] Yuji Nakatsukasa and Joel A. Tropp.  Fast & accurate randomized algorithms for linear systems and eigenvalue problems.  arXiv preprint arXiv:2111.00113, 2021.
  * [25] Balas Kausik Natarajan.  Sparse approximate solutions to linear systems.  SIAM J. Comput., 24(2):227–234, 1995.
  * [26] Yagyensh Chandra Pati, Ramin Rezaiifar, and Perinkulam Sambamurthy Krishnaprasad.  Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition.  In Proceedings of 27th Asilomar Conference on Signals, Systems and Computers, pages 40–44. IEEE, 1993.
  * [27] Bernard Philippe and Lothar Reichel.  On the generation of Krylov subspace bases.  Appl. Numer. Math., 62(9):1171–1186, 2012.
  * [28] Yousef Saad.  Variations on Arnoldi’s method for computing eigenelements of large unsymmetric matrices.  Linear Algebra Appl., 34:269–295, 1980.
  * [29] Yousef Saad.  Krylov subspace methods for solving large unsymmetric linear systems.  Math. Comput., 37(155):105–126, 1981.
  * [30] Yousef Saad.  Iterative Methods for Sparse Linear Systems, 2nd edition.  SIAM, Philadelphia, 2000.
  * [31] Yousef Saad.  Numerical Methods for Large Eigenvalue Problems: Revised Edition.  SIAM, 2011.
  * [32] Tamas Sarlos.  Improved approximation algorithms for large matrices via random projections.  In 47th Annual IEEE Symposium on Foundations of Computer Science (FOCS’06), pages 143–152. IEEE, 2006.
  * [33] Henrik V Sorensen and C Sidney Burrus.  Efficient computation of the DFT with only a subset of input or output points.  IEEE Trans. Signal Process., 41(3):1184–1200, 1993.
  * [34] Andreas Stathopoulos, Yousef Saad, and Kesheng Wu.  Dynamic thick restarting of the Davidson and the implicitly restarted Arnoldi methods.  SIAM J. Sci. Comput., 19(1):227–245, 1998.
  * [35] G. W. Stewart.  A Krylov–Schur algorithm for large eigenproblems.  SIAM J. Matrix Anal. Appl., 23(3):601–614, 2001/02.
  * [36] Robert Tibshirani.  The LASSO method for variable selection in the Cox model.  Stat. Med., 16(4):385–395, 1997.
  * [37] Joel A. Tropp.  Greed is good: Algorithmic results for sparse approximation.  IEEE Trans. Inform. Theory, 50(10):2231–2242, 2004.
  * [38] David P. Woodruff.  Sketching as a tool for numerical linear algebra.  Found. Trends Theor. Comput. Sci., 10(1–2):1–157, 2014.
  * [39] Franco Woolfe, Edo Liberty, Vladimir Rokhlin, and Mark Tygert.  A fast randomized algorithm for the approximation of matrices.  Appl. Comput. Harmon. Anal., 25:335–366, 2008.