remarkRemark hypothesisHypothesis claimClaim

# GMRES with randomized sketching and deflated restarting

Liam Burke School of Mathematics, Trinity College Dublin, College Green,
Dublin 2, Ireland, . L. B. acknowledges funding from an Irish Research Council
Government of Ireland Postgraduate Scholarship<EMAIL_ADDRESS>Stefan Güttel.
S. G.’s work is partly funded by a Royal Society Industry Fellowship
IF/R1/231032 Department of Mathematics, The University of Manchester, M13 9PL
Manchester, United Kingdom<EMAIL_ADDRESS>Kirk M. Soodhalter
School of Mathematics, Trinity College Dublin, College Green, Dublin 2,
Ireland<EMAIL_ADDRESS>
###### Abstract

We present a new Krylov subspace recycling method for solving a linear system
of equations, or a sequence of slowly changing linear systems. Our new method,
named GMRES-SDR, combines randomized sketching and deflated restarting in a
way that avoids orthogononalizing a full Krylov basis. We provide new theory
which characterizes sketched GMRES with and without augmentation as a
projection method using a semi-inner product. We present results of numerical
experiments demonstrating the effectiveness of GMRES-SDR over competitor
methods such as GMRES-DR and GCRO-DR.

###### keywords:

linear systems, Krylov method, subspace recycling, randomized sketching

65F60, 65F50, 65F10, 68W20

##  1 Introduction

This paper is concerned with the development of a Krylov subspace recycling
algorithm for the efficient solution of a linear system of equations

$A\mathbf{x}=\mathbf{b}$

with a nonsymmetric matrix $A\in\mathbb{C}^{N\times N}$ and a right-hand side
vector $\mathbf{b}\in\mathbb{C}^{N}$. Our algorithm can also naturally be
applied to solve sequences of nearby linear systems

(1) $A^{(i)}\mathbf{x}^{(i)}=\mathbf{b}^{(i)},\qquad i=1,2,\ldots,$

using accumulated Krylov information computed to improve convergence as the
sequence progresses. Such problems arise in a variety of scientific computing
applications (see, e.g, [20]), and Krylov subspace recycling has been
demonstrated to be a successful technique in minimizing the computational
resources and runtime required to solve the full sequence of problems [5, 12].
Recycling algorithms belong to the class of augmented Krylov subspace methods,
where the augmentation subspace for each problem is constructed or recycled
from the Krylov subspace used to solve a previous problem in the sequence.
Although recycling introduces additional computational overhead per iteration
when compared to non-augmented methods, the presence of the recycling subspace
can often allow each problem to converge faster and in much less iterations,
resulting in a reduction in overall computational cost and runtime for the
full problem sequence.

One of the dominating computational kernels of common restarting or recycling
methods, such as GMRES-DR [15] or GCRO-DR [20], is the orthogonalization of a
full Krylov basis in each cycle. In this paper, we build on existing work in
[16, 2, 13, 7] to develop a new GMRES-based recycling method, named GMRES-SDR,
which reduces the expensive orthogonalization costs using the well established
technique of randomized sketching [26].

The structure of the paper is as follows. In Section 2 we briefly review some
of the most commonly used recycling methods and introduce our notation.
Section 3 introduces the new GMRES-SDR method, including the implementation
details. As analysis of augmented and recycled Krylov subspace methods is
often built on top of the analysis of the method being augmented, in Section 4
we develop some additional theory about sketched GMRES as a projection-type
scheme and its behavior. We build on this in Section 5 by showing that
sketched augmented GMRES exhibits many of the features of an augmented
projection method as laid out in [25], and we extend the analysis of sketched
GMRES developed in Section 4 to get some computable convergence bounds for
sketched augmented GMRES. In Section 6, we demonstrate that the sketched
augmented GMRES method we propose in this manuscript is effective and reduces
time-to-solution over other augmented and recycled methods for some large-
scale problems arising in the computational sciences.

##  2 Relevant background

This paper concerns augmented and sketched varieties of GMRES, and so we need
to introduce some relevant theory and notation concerning GMRES and
augmented/recycled varieties thereof. This discussion is abbreviated, focusing
on the aspects we need to develop and analyze the methods we propose in
subsequent sections.

###  2.1 The classic GMRES method

We describe one cycle of classic GMRES [23] for solving a single linear system
$A\mathbf{x}=\mathbf{b}$ with initial guess $\mathbf{x}_{0}$ and associated
residual $\mathbf{r}_{0}:=\mathbf{b}-A\mathbf{x}_{0}$. GMRES generates an
orthonormal basis for the Krylov space
$\mathcal{K}_{m}(A,\mathbf{r}_{0}):=\operatorname{span}\\{\mathbf{r}_{0},A\mathbf{r}_{0},\ldots,A^{m-1}\mathbf{r}_{0}\\}$,
which are stored as the columns of a matrix $V_{m}\in\mathbb{C}^{N\times m}$.
This basis is generated by the modified Gram–Schmidt process, satisfying the
Arnoldi relation

(2) $\displaystyle AV_{m}=V_{m+1}\underline{H_{m}},$

where $\underline{H_{m}}\in\mathbb{C}^{(m+1)\times m}$ is upper Hessenberg.
GMRES computes the minimizer

$\displaystyle\mathbf{x}_{m}=\underset{\mathbf{t}\in\mathcal{K}_{m}(A,\mathbf{r}_{0})}{\mathrm{argmin}}\left\|\mathbf{r}_{0}-A\mathbf{t}\right\|,$

which is equivalent to setting $\mathbf{x}_{m}=V_{m}\mathbf{y}_{m}$ with

$\displaystyle\mathbf{y}_{m}=\underset{\mathbf{y}\in\mathbb{C}^{m}}{\mathrm{argmin}}\left\|\underline{H_{m}}\mathbf{y}-\left\|\mathbf{r}_{0}\right\|\mathbf{e}_{1}\right\|.$

The GMRES least squares approximation for $A\mathbf{x}=\mathbf{b}$ using any
basis $V_{m}\in\mathbb{C}^{N\times m}$ can be written as
$\mathbf{x}_{m}=\mathbf{x}_{0}+V_{m}\mathbf{y}_{m}$, where
$\mathbf{y}_{m}=\left(V_{m}^{\ast}A^{\ast}AV_{m}\right)^{-1}V_{m}^{\ast}A^{\ast}\mathbf{r}_{0}$.
The residual can thus be written as

$\displaystyle\mathbf{r}_{m}=(I-AV_{m}\left(V_{m}^{\ast}A^{\ast}AV_{m}\right)^{-1}V_{m}^{\ast}A^{\ast})\mathbf{r}_{0},$

which is an orthogonal projection of $\mathbf{r}_{0}$ onto
$\mathrm{Range}(AV_{m})^{\perp}$, leading to the well-known fact that
$\mathrm{Range}(AV_{m})$ is the GMRES residual constraint space. Let
$\Phi_{\mathcal{K}_{m}}=AV_{m}\left(V_{m}^{\ast}A^{\ast}AV_{m}\right)^{-1}V_{m}^{\ast}A^{\ast}$
be the orthogonal projector onto $\mathrm{Range}(AV_{m})$, which we call a
_residual projector_. There is a sibling _error projector,_
$\Pi_{\mathcal{K}_{m}}=V_{m}\left(V_{m}^{\ast}A^{\ast}AV_{m}\right)^{-1}V_{m}^{\ast}A^{\ast}A$,
which is the $A^{\ast}A$-orthogonal projector onto the space
$\mathcal{K}_{m}(A,\mathbf{r}_{0})$ so that
$\mathbf{x}_{m}=\mathbf{x}_{0}+\Pi_{\mathcal{K}_{m}}\boldsymbol{\eta}_{0}$,
where $\boldsymbol{\eta}_{0}=\mathbf{x}-\mathbf{x}_{0}$.

###  2.2 Recycling and augmentation of minimum residual methods

In [25] it was established that augmented subspace methods (such as augmented
Krylov methods and recycling methods) can all be characterized in a general
_projection–correction_ framework. In the case of a minimum residual method,
choosing $\mathbf{f}_{m}\in\mathrm{Range}\left(\mathbf{V}_{m}\right)$ for
$\mathbf{V}_{m}=[U,V_{m}]$ such that

$\displaystyle\mathbf{f}_{m}=\underset{\mathbf{f}\in\mathrm{Range}\left(\mathbf{V}_{m}\right)}{\mathrm{argmin}}\left\|\mathbf{b}-A\left(\mathbf{x}_{0}+\mathbf{f}\right)\right\|$

is equivalent to selecting

$\displaystyle\mathbf{t}_{m}=\underset{\mathbf{t}\in\mathrm{Range}\left(V_{m}\right)}{\mathrm{argmin}}\left\|\left(I-\Phi_{U}\right)\left(\mathbf{b}-A\mathbf{t}\right)\right\|,$

where $\Phi_{U}$ is the orthogonal residual projector onto
$\mathrm{Range}(AU)$, and constructing
$\mathbf{x}_{m}=\mathbf{x}_{0}+\Pi_{U}\boldsymbol{\eta}_{0}+(I-\Pi_{U})\mathbf{t}_{m}$,
where $\Pi_{U}$ is the sibling $A^{\ast}A$-orthogonal error projector onto
$\mathrm{Range}\left(U\right)$. In other words, the augmented minimum residual
method approach is equivalent to applying a minimum residual method to the
projected subproblem

(3)
$\displaystyle\left(I-\Phi_{U}\right)A\mathbf{t}=\left(I-\Phi_{U}\right)\mathbf{b}$

over the approximation space $\mathrm{Range}\left(V_{m}\right)$ and then
applying some projections back into $\mathrm{Range}\left(U\right)$ to
construct $\mathbf{x}_{m}$.

###  2.3 GMRES-DR and GCRO-DR are sometimes equivalent

It was observed in [19] that GCRO-DR [20] and GMRES-DR [15] are mathematically
equivalent in limited circumstances.

The GCRO-DR method [20] takes
$\mathrm{Range}\left(V_{m}\right)=\mathcal{K}_{m}\left(\left(I-\Phi_{U}\right)A,\left(I-\Phi_{U}\right)\mathbf{r}_{0}\right)$
with $V_{m}$ being built by the Arnoldi process. Thus it satisfies the Arnoldi
relation for the projected operator,

$\displaystyle\left(I-\Phi_{U}\right)AV_{m}=V_{m+1}\underline{H_{m}}.$

From this, it follows that the minimization is equivalent to applying a GMRES
iteration to (3). This method accommodates beginning the iteration with an
arbitrary $U$, and any technique of selecting vectors to recycle can be used.

The GMRES-DR method [15] takes a similar approach but is restricted to reusing
subspace information for a single linear system between restarts. At the end
of a cycle, the method computes harmonic Ritz vectors which are then
orthogonalized to save for the next cycle. Let $U$ denote the matrix with the
harmonic Ritz vectors as columns. The last GMRES residual $\mathbf{r}_{m}$ is
orthogonalized to produce $\mathbf{v}_{k+1}=(I-UU^{\ast})\mathbf{r}_{m}$.
These orthogonalizations are performed efficiently by taking advantage of the
Arnoldi relations of the just-completed GMRES cycle. It is proven that if
harmonic Ritz vectors are used, the space
$\mathrm{span}\left\\{\mathbf{u}_{1},\mathbf{u}_{2},\ldots,\mathbf{u}_{k},\mathbf{v}_{k+1}\right\\}$
is a Krylov subspace and an Arnoldi relation for it can be cheaply built.
Thereafter, the iteration continues. Because the augmentation space
$\mathrm{Range}\left(U\right)$ can be generated so that the augmented subspace
remains a Krylov subspace, GMRES-DR has a leaner implementation than GCRO-DR.
In particular, with a proper Arnoldi relation, only $U$ needs to be stored.
There is no need to store $AU$.

Consider running a cycle of GMRES, stopping to restart, and computing some
harmonic Ritz vectors to carry over to the next cycle. GMRES-DR proceeds in
the next cycle by constructing an orthonormal basis for $\mathbf{V}_{m}$ which
is established to be a Krylov subspace (since $U$ consists of harmonic Ritz
vectors from the previous cycle). We observe that $\mathbf{v}_{j}$ for
$j=k+2,k+3,\cdots,m$ is generated by applying $A$ to $\mathbf{v}_{j-1}$,
orthogonalizing against the columns of $U$, and orthogonalizing against
$\mathbf{v}_{i}$ for $i=k+1,k+2,\ldots,j-1$, i.e., perform the Arnoldi
process. This means that
$\left\\{\mathbf{v}_{k+1},\mathbf{v}_{k+2},\ldots,\mathbf{v}_{j-1},\mathbf{v}_{j}\right\\}$
forms a basis for
$\mathcal{K}_{j-k-1}\left((I-UU^{\ast})A,(I-UU^{\ast})\mathbf{r}_{m}\right)$,
the space which would be built by GCRO-DR in this setting. One observes this
by noting that $U$ being built from harmonic Ritz vectors means they satisfy
an Arnoldi relation $AU=[U,\mathbf{v}_{k+1}]\underline{\widetilde{H}_{m}},$
where Morgan [15] derives a cheap procedure for obtaining the upper
Hesssenberg matrix $\underline{\widetilde{H}_{m}}$.

##  3 Derivation of GMRES-SDR

Let us first consider a single linear system $A\mathbf{x}=\mathbf{b}$ with
initial guess $\mathbf{x}_{0}$ and associated residual
$\mathbf{r}_{0}:=\mathbf{b}-A\mathbf{x}_{0}$. Assume that
$\mathbf{V}_{m}=[U,V_{m}]$ spans some search space, typically with the columns
of $U$ approximating an invariant subspace for the matrix $A$, and $V_{m}$
being a (not necessarily orthogonal) Krylov basis of
$\mathcal{K}_{m}(A,\mathbf{r}_{0})=\operatorname{span}\\{\mathbf{r}_{0},A\mathbf{r}_{0},\ldots,A^{m-1}\mathbf{r}_{0}\\}$.

In order to exploit randomized sketching (see, e.g., [14, 3, 4, 1, 16]), we
further assume that we have an operator $S\in\mathbb{C}^{s\times N}$ with
$m<s\ll N$ which acts as an approximate isometry for the Euclidean norm
$\|\,\cdot\,\|$. More precisely, given a positive integer $m$ and some
$\varepsilon\in[0,1)$, let $S$ be such that for all vectors
$\mathbf{v}\in\operatorname{span}(\mathbf{V}_{m})$,

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

The mapping $S$ is called an _$\varepsilon$ -subspace embedding_ for
$\operatorname{span}(\mathbf{V}_{m})$; see, e.g., [24, 26, 14]. Condition (4)
can equivalently be stated with the Euclidean inner product [24, Cor. 4]: for
all $\mathbf{u},\mathbf{v}\in\operatorname{span}(\mathbf{V}_{m})$,

${}\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, $S$ is not explicitly available but we can draw it at random to
achieve (4) with high probability.

Then by [7, Remark 3.1], the augmented GMRES-type approximant for a matrix
function $f(A)\mathbf{b}$ is

(5)
$\widetilde{\mathbf{x}}_{m}=\mathbf{V}_{m}f\left([(S\mathbf{W}_{m})^{*}S\mathbf{V}_{m}]^{-1}(S\mathbf{W}_{m})^{*}S\mathbf{W}_{m}\right)[(S\mathbf{W}_{m})^{*}S\mathbf{V}_{m}]^{-1}(S\mathbf{W}_{m})^{*}S\mathbf{r}_{0},$

with $\mathbf{W}_{m}=A\mathbf{V}_{m}$. Here we have $f(z)=z^{-1}$, in which
case this simplifies to

$\widetilde{\mathbf{x}}_{m}=\mathbf{V}_{m}\left[(SA\mathbf{V}_{m})^{*}SA\mathbf{V}_{m}\right]^{-1}(SA\mathbf{V}_{m})^{*}S\mathbf{r}_{0},$

or alternatively,

(6)
$\widetilde{\mathbf{x}}_{m}=\mathbf{V}_{m}\widetilde{\mathbf{y}}_{m},\quad\widetilde{\mathbf{y}}_{m}=(SA\mathbf{V}_{m})^{\dagger}(S\mathbf{r}_{0}).$

This means that

$\widetilde{\mathbf{x}}_{m}\ \text{minimizes}\
\texttt{sres}:=\|S\mathbf{r}_{0}-SA\widetilde{\mathbf{x}}\|\ \text{over all}\
\widetilde{\mathbf{x}}\in\operatorname{span}(\mathbf{V}_{m}).$

Note that this is also the definition of the sGMRES approximant proposed in
[16, eq. (1.6)]. The interesting aspect here is that this approximant does not
rely in orthogonality of $\mathbf{V}_{m}$, yet its residual
$\widetilde{\mathbf{r}}_{m}=\mathbf{r}_{0}-A\widetilde{\mathbf{x}}_{m}$ is
closely related to the residual
$\mathbf{r}_{m}=\mathbf{r}_{0}-A\mathbf{x}_{m}$ of the full GMRES approximant
$\mathbf{x}_{m}:=\mathbf{V}_{m}(A\mathbf{V}_{m})^{\dagger}\mathbf{r}_{0}$.
This latter approximant

${\mathbf{x}}_{m}\ \text{minimizes}\
\|\mathbf{r}_{0}-A\widetilde{\mathbf{x}}\|\ \text{over all}\
\widetilde{\mathbf{x}}\in\operatorname{span}(\mathbf{V}_{m}).$

Combining this with (4), we have the following residual relations

(7)
$\|\mathbf{r}_{m}\|\leq{\|\widetilde{\mathbf{r}}_{m}\|}\leq\frac{1}{\sqrt{1-\varepsilon}}\|S\widetilde{\mathbf{r}}_{m}\|\leq\frac{1}{\sqrt{1-\varepsilon}}\|S\mathbf{r}_{m}\|\leq\sqrt{\frac{1+\varepsilon}{1-\varepsilon}}\|\mathbf{r}_{m}\|.$

###  3.1 Implementation

The GMRES-SDR method is now fairly straightforward to implement, as detailed
in Algorithm 1. There are several components that differ from the usual GMRES
implementation and we discuss them separately in the following subsections.

Algorithm 1 One cycle of GMRES-SDR for $A\mathbf{x}=\mathbf{r}_{0}$

1: Input: ${A}\in\mathbb{C}^{N\times N}$, $\mathbf{r}_{0}\in\mathbb{C}^{N}$,
$U\in\mathbb{C}^{N\times\widehat{k}}$, $S\in\mathbb{C}^{s\times N}$, tol,
integers $m$, $t$, $k$

2: Optional input: $SU\in\mathbb{C}^{s\times\widehat{k}}$,
$SAU\in\mathbb{C}^{s\times\widehat{k}}$, $\texttt{safety}\geq 1$ (default 1.4)

3: Output: Approximant $\widetilde{\mathbf{x}}\approx A^{-1}\mathbf{r}_{0}$
and recycling subspace $U,SU,SAU$

4: If not provided as input, compute sketches $SU$ and $SAU$

5: Sketch initial residual $S\mathbf{r}_{0}$

6: $V:=[\mathbf{r}_{0}/\|\mathbf{r}_{0}\|,O_{N\times m}]$, $H:=O_{(m+1)\times
m}$

7: $SV:=[S\mathbf{r}_{0}/\|\mathbf{r}_{0}\|,O_{N\times m}]$, $SAV:=O_{s\times
m}$

8: % —— Truncated Arnoldi process with QR-based residual estimation ——

9: for $j=1,2,\ldots,m$

10: $\mathbf{w}:=AV(\,:\,,j)$

11: for $i=\max(j-t+1,1)$

12: $H(i,j):=V(\,:\,,i)^{*}\mathbf{w}$

13: $\mathbf{w}:=\mathbf{w}-V(\,:\,,i)H(i,j)$

14: $H(j+1,j):=\|\mathbf{w}\|$

15: $V(\,:\,,j+1):=\mathbf{w}/H(j+1,j)$

16: $SV(\,:\,,j+1):=S\cdot V(\,:\,,j+1)$ % one sketch per iter

17: $SAV(\,:\,,j):=SV(\,:\,,1:j+1)H(1:j+1,j)$

18: $SAW:=[SAU,SAV(\,:\,,1:j)]$

19: Compute (or update) economic $[Q,R]=\texttt{qr}(SAW)$

20: Compute least squares solution
$\widetilde{\mathbf{y}}:=R^{-1}Q^{*}[S\mathbf{r}_{0}]$

21: $\texttt{sres}:=\|S\mathbf{r}_{0}-SAW\cdot\widetilde{\mathbf{y}}\|$ %
sketched residual

22: if $\texttt{sres}<\texttt{tol}/\texttt{safety}$ or $j=m$

23: $\widetilde{\mathbf{x}}:=[U,V(\,:\,,1:j)]\cdot\widetilde{\mathbf{y}}$

24: $\texttt{res}:=\|\mathbf{r}_{0}-A\widetilde{\mathbf{x}}\|$ % true residual

25: if $\texttt{res}<\texttt{tol}$

26: break % stop trunc. Arnoldi

27: else

28: $\texttt{safety}:=\texttt{res}/\texttt{sres}$ % increase safety factor

29: % ——— SVD-based update of recycling space using harmonic Ritz ———

30: Compute truncated economic SVD
$\widehat{U}_{\ell}\Sigma_{\ell}\widehat{V}_{\ell}^{*}\approx SAW$

31: Compute $M_{\ell}:=\widehat{U}_{\ell}^{*}\cdot SW\cdot\widehat{V}_{\ell}$

32: Compute ordered QZ decomposition
$[Q,Z]=\texttt{qz}(M_{\ell},\Sigma_{\ell})$

33: Compute $U:=[U,V(\,:\,,1:m)]\widehat{V}_{\ell}Z(\,:\,,1:k)$

34: Compute $SU:=[SU,SV(\,:\,,1:m)]\widehat{V}_{\ell}Z(\,:\,,1:k)$

35: Compute $SAU:=[SAU,SAV(\,:\,,1:m)]\widehat{V}_{\ell}Z(\,:\,,1:k)$

36: return $\widetilde{\mathbf{x}}$, $U$, $SU$, $SAU$

###  3.2 Krylov basis generation

In order to generate the Krylov bases $V_{m}$ at a cost that grows only
linearly with $m$, we use the truncated Arnoldi process in Algorithm 1. This
means that each new basis vector is projected against the previous $t\ll m$
vectors only.

###  3.3 Solution of the projected least squares problem

In standard GMRES the projected least squares problem is usually solved by
updating a QR factorization of the Hessenberg matrix $\underline{H_{m}}$ using
Givens rotations. Here, the matrix $SA\mathbf{V}_{m}$ in the least squares
problem (6) is not of Hessenberg form. However, the matrix
$\mathbf{V}_{m}=[U,V_{m}]$ grows one column at a time when the truncated
Arnoldi process for generating $V_{m}$ progresses, hence we can still update
the QR factorization of $SA\mathbf{V}_{m}$ when a column is added, e.g., by
using Householder-QR or the modified Gram–Schmidt with reorthogonalization.

###  3.4 Residual control and updating safety

The residual of the sketched problem (6) can serve as an approximation for the
true residual of the current approximant using

${\|\widetilde{\mathbf{r}}_{m}\|}\leq\frac{1}{\sqrt{1-\varepsilon}}\|S\widetilde{\mathbf{r}}_{m}\|.$

However, we generally do not know what $\varepsilon$ is. Hence we propose the
following strategy: within the GMRES-SDR cycle we monitor the sketched
residual norm $\|S\widetilde{\mathbf{r}}_{m}\|$. Only if this norm satisfies

$\|S\widetilde{\mathbf{r}}_{m}\|<\texttt{tol}/\texttt{safety}$

for some safety factor $\texttt{safety}>1$, we compute the true residual norm
$\|\widetilde{\mathbf{r}}_{m}\|$. If also this quantify is smaller than tol,
we can terminate the Arnoldi iteration. Otherwise we will increase the safety
factor to

$\texttt{safety}=\|\widetilde{\mathbf{r}}_{m}\|/\|S\widetilde{\mathbf{r}}_{m}\|$

and continue.

An example illustration is given in Figure 1. The restart length in this
example is chosen as $m=100$, and it is interesting to note that the
distortion $\|\widetilde{\mathbf{r}}_{m}\|/\|S\widetilde{\mathbf{r}}_{m}\|$
follows a generally increasing trend within each cycle. We also notice a spike
in the measured distortion of the basis vectors,
$\|\mathbf{v}_{j}\|/\|S\mathbf{v}_{j}\|=1/\|S\mathbf{v}_{j}\|$, with the
height of the spike being very close to the distortion on residual. This is
likely caused by the GMRES property that the Krylov basis vector
$\mathbf{v}_{1}$ at the beginning of a cycle is strongly correlated with the
last basis vector of the previous cycle (in FOM these vectors would be
perfectly collinear). In our experiments we found that an initial value
$\texttt{safety}=1.4$ works well, with the computation of the true residual
only needed very rarely (apart from the unavoidable computation at the end of
each Arnoldi cycle).

Figure 1: Example showing measured sketching distortion factors as the GMRES-
SDR iteration progresses, as well as the ratio between the true and sketched
residuals. In this example, the Krylov dimension of size $m=100$, a truncation
parameter of $t=2$, a recycling subspace dimension of size $k=20$, and a
sketching parameter $s=500$ was used.

###  3.5 Harmonic updating of the augmentation space

We now turn our attention to the problem of updating the augmentation space
$U$ from one restart cycle or problem to the next. Say, in one cycle we have
computed a Krylov basis $V_{m}$ of $\mathcal{K}_{m}(A,\mathbf{r}_{0})$ and we
have an augmentation basis $U\in\mathbb{C}^{N\times k}$ accumulated over
previous cycles. Neither $V_{m}$ nor $U$ are now assumed to be orthonormal.

Consider the augmented basis matrix $\mathbf{V}_{m}=[U,V_{m}]$ and assume that
we have the sketches $S\mathbf{V}_{m}$ and $SA\mathbf{V}_{m}$ at our disposal.
The (non-sketched) harmonic Ritz pairs [17] are defined as
$(\vartheta,\mathbf{z}=\mathbf{V}_{m}\mathbf{y}_{m})$ such that

$A\mathbf{z}-\vartheta\mathbf{z}\perp A\mathbf{V}_{m}.$

Therefore,

$(A\mathbf{V}_{m})^{*}(A\mathbf{V}_{m}\mathbf{y}_{m}-\vartheta\mathbf{V}_{m}\mathbf{y}_{m})=0,$

which means that

$(A\mathbf{V}_{m})^{*}(A\mathbf{V}_{m})\mathbf{y}_{m}=\vartheta(A\mathbf{V}_{m})^{*}\mathbf{V}_{m}\mathbf{y}_{m},$

or equivalently,

$\mathbf{y}_{m}=\vartheta[(A\mathbf{V}_{m})^{*}(A\mathbf{V}_{m})]^{-1}(A\mathbf{V}_{m})^{*}\mathbf{V}_{m}\mathbf{y}_{m}=\vartheta(A\mathbf{V}_{m})^{\dagger}\mathbf{V}_{m}\mathbf{y}_{m}.$

In other words, the harmonic Ritz pairs can be computed by solving an
eigenvalue problem

$M\mathbf{y}_{m}:=[(A\mathbf{V}_{m})^{\dagger}\mathbf{V}_{m}]\mathbf{y}_{m}=\vartheta^{-1}\mathbf{y}_{m}.$

Now to utilize sketching, we realize that each column of the matrix $M$ can be
thought of as the solution of a least squares problem with the tall skinny
matrix $A\mathbf{V}_{m}$. The sketched counterpart of this is

$\widehat{M}\widehat{\mathbf{y}}_{m}:=[(SA\mathbf{V}_{m})^{\dagger}S\mathbf{V}_{m}]\widehat{\mathbf{y}}_{m}=\widehat{\vartheta}^{-1}\widehat{\mathbf{y}}_{m}.$

Note that $\widehat{M}$ is the solution (of minimum Frobenius norm) to

(8)
$\widehat{M}:=\arg\min_{M\in\mathbb{C}^{(m+k)\times(m+k)}}\|SA\mathbf{V}_{m}M-S\mathbf{V}_{m}\|_{F},$

a harmonic version of the sRR approach presented in [16, Sec. 6.3].

To stabilize the harmonic Ritz extraction, we exploit a truncated SVD
$SA\mathbf{V}_{m}\approx\widehat{U}_{\ell}\Sigma_{\ell}\widehat{V}_{\ell}^{*}$
and replace

$SA\mathbf{V}_{m}\rightarrow\widehat{U}_{\ell}\quad\text{and}\quad
S\mathbf{V}_{m}\rightarrow(S\mathbf{V}_{m})\widehat{V}_{\ell}\Sigma_{\ell}^{-1},$

leading to the regularized

$\widehat{M}_{\ell}:=\widehat{U}_{\ell}^{*}(S\mathbf{V}_{m})\widehat{V}_{\ell}\Sigma_{\ell}^{-1}\approx\widehat{M}.$

We then compute an ordered QZ factorization such that

$Q[\widehat{U}_{\ell}^{*}(S\mathbf{V}_{m})\widehat{V}_{\ell}]Z=Q[\Sigma_{\ell}]Z$

is upper triangular with the $k$ most relevant generalized eigenvalues
appearing in the top-left block. As we are interested in harmonic Ritz values,
we should target the _largest modulus_ eigenvalues. The desired harmonic Ritz
vectors are then given as

$U:=[U,V_{m}]\widehat{V}_{\ell}Z(\,:\,,1:k).$

Similarly, the sketched version of $U$ and $AU$ can be computed without
explicit sketches as shown in Algorithm 1.

###  3.6 Solving slowly varying linear systems

So far, we have introduced GMRES-SDR for a single linear system
$A\mathbf{x}=\mathbf{b}$. Let us now consider the cases that
$A\mathbf{x}^{(i)}=\mathbf{b}^{(i)}$ with varying right-hand side vectors, and
$A^{(i)}\mathbf{x}^{(i)}=\mathbf{b}^{(i)}$ with varying vectors and system
matrices.

In the first case, when only $\mathbf{b}^{(i)}$ varies, we can simply reuse
the recycling subspace $U,SU,SAU$ from one problem to the next. We just take
the outputs $U,SU,SAU$ from Algorithm 1 computed during the solution of
$A\mathbf{x}^{(i)}=\mathbf{b}^{(i)}$ and feed them as inputs to Algorithm 1
when solving $A\mathbf{x}^{(i+1)}=\mathbf{b}^{(i+1)}$.

The second case with problem-dependent system matrices,
$A^{(i)}\mathbf{x}^{(i)}=\mathbf{b}^{(i)}$, allows for two variants which we
refer to _exact_ and _inexact_ GMRES-SDR, respectively. In the exact variant
we reuse $U,SU$ from one problem $i$ to problem $i+1$, but compute
$SA^{(i+1)}U$ by explicit matrix multiplication and sketching. This requires
$k$ additional matrix-vector products and the sketching of $k$ additional
vectors for each new problem.

In the inexact GMRES-SDR variant we feed $SA^{(i)}U$ as an input to Algorithm
1 even though we are now solving
$A^{(i+1)}\mathbf{x}^{(i+1)}=\mathbf{b}^{(i+1)}$ with a different system
matrix. As a consequence, the matrix $SAW$ computed in line 18 of Algorithm 1
becomes

$SAW=[SA^{(i)}U,SA^{(i+1)}V(:,1:j)],$

i.e., we are mixing Krylov spaces of two different matrices. In particular,
the least squares problem in line 21 no longer corresponds to solving a
sketched version of the full GMRES problem for
$A^{(i+1)}\mathbf{x}^{(i+1)}=\mathbf{b}^{(i+1)}$ and the residual relations
(7) are violated. Still, we have observed numerically that one can sometimes
get away with the inexactness when $A^{(i)}$ changes only very slightly from
one problem to the next. However, there is no guarantee that this works and
the sketched residuals may even increase as the iteration progresses. We will
demonstrate this in Section 6.3.

##  4 Analysis of sketched GMRES

We build on the one-cycle description of sketched GMRES as introduced in [16].
They describe a general matrix $V_{m}\in\mathbb{C}^{N\times m}$ whose columns
are some basis vectors of a Krylov subspace, i.e.,
$\mathrm{Range}(V_{m})=\mathcal{K}_{m}(A,\mathbf{r}_{0})$. We extend existing
GMRES convergence theory to the sketched GMRES setting, separating into basis-
agnostic and basis-dependent results. Then in Section 5 we put sketched
augmented GMRES into the existing framework of augmented iterative methods.
This allows us to extend the analysis for sketched GMRES to augmented sketched
GMRES.

###  4.1 Sketching the constraint space for one cycle of sketched GMRES

We pick up from the discussion in Section 2.1 and take a similar approach to
understand the constraint for sketched GMRES. For the same basis $V_{m}$, the
sketched GMRES approximation is defined as

$\displaystyle\widetilde{\mathbf{x}}_{m}=\mathbf{x}_{0}+V_{m}\widetilde{\mathbf{y}}_{m},$

where $\widetilde{\mathbf{y}}_{m}$ minimizes
$\left\|S(AV_{m}\mathbf{y}-\mathbf{r}_{0})\right\|$, and via the normal
equations, it follows that

$\displaystyle\widetilde{\mathbf{y}}_{m}=\left(V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}\right)^{\dagger}\left(V_{m}^{\ast}A^{\ast}S^{\ast}S\right)\mathbf{r}_{0}.$

(We establish in Lemma 4.3 below that $V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}$ is
invertible for a sufficiently small sketching parameter $\varepsilon$.) The
sketched GMRES residual satisfies the relation

$\displaystyle\widetilde{\mathbf{r}}_{m}=\mathbf{r}_{0}-AV_{m}\left(V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}\right)^{\dagger}\left(V_{m}^{\ast}A^{\ast}S^{\ast}S\right)\mathbf{r}_{0}=\left(I-\widehat{\Phi}_{\mathcal{K}_{m}}\right)\mathbf{r}_{0},$

where
$\widehat{\Phi}_{\mathcal{K}_{m}}=AV_{m}\left(V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}\right)^{\dagger}\left(V_{m}^{\ast}A^{\ast}S^{\ast}S\right)$.
Similarly, if we let

$\widehat{\Pi}_{\mathcal{K}_{m}}=V_{m}\left(V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}\right)^{\dagger}\left(V_{m}^{\ast}A^{\ast}S^{\ast}S\right)A,$

the sketched GMRES approximation satisfies
$\hat{\mathbf{x}}_{m}=\mathbf{x}_{0}+\widehat{\Pi}_{\mathcal{K}_{m}}\boldsymbol{\eta}_{0}$.
We observe that $S\widehat{\Phi}_{\mathcal{K}_{m}}=\Phi_{SAV_{m}}S$, where
$\Phi_{SAV_{m}}$ is the orthogonal projector onto
$\mathrm{Range}\left(SAV_{m}\right)$, meaning that the sketched residual
satisfies $S\hat{\mathbf{r}}_{m}\perp\mathrm{Range}\left(SAV_{m}\right)$. This
is easy to see since we solved a sketched least squares problem.

###  4.2 Sketched GMRES is a projection method

The matrices $\widehat{\Pi}_{\mathcal{K}_{m}}$ and
$\widehat{\Phi}_{\mathcal{K}_{m}}$ serve the same role as the error and
residual projectors in GMRES. We observe that they actually have many
properties of projectors. This builds on ideas in, e.g., [3, 18]. Note that in
the following discussion, the space spanned by the columns of
$S^{\ast}SAV_{m}$ is discussed. Since $S$ sketches the true GMRES constraint
space, the space must be lifted or interpolated back up to the full space in
order to be used as a sketched residual constraint. The operator $S^{\ast}$
lifts the sketched space back to a subspace in $n$ dimensions.

######  Lemma 4.1.

The matrices $\widehat{\Pi}_{\mathcal{K}_{m}}$ and
$\widehat{\Phi}_{\mathcal{K}_{m}}$ are both idempotent, with
$\mathrm{Range}\left(\widehat{\Phi}_{\mathcal{K}_{m}}\right)\subseteq\mathrm{Range}(AV_{m})$
with null space containing $\mathrm{Range}(S^{\ast}SAV_{m})^{\perp}$ and
$\mathrm{Range}\left(\widehat{\Pi}_{\mathcal{K}_{m}}\right)\subseteq\mathrm{Range}(V_{m})$
with null space containing $\mathrm{Range}(A^{\ast}S^{\ast}SAV_{m})^{\perp}$.

######  Proof 4.2.

Idemptoency can be proven directly. For $\widehat{\Phi}_{\mathcal{K}_{m}}$, we
can write

$\displaystyle\widehat{\Phi}_{\mathcal{K}_{m}}^{2}$
$\displaystyle=AV_{m}\left(V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}\right)^{\dagger}\left(V_{m}^{\ast}A^{\ast}S^{\ast}S\right)AV_{m}\left(V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}\right)^{\dagger}\left(V_{m}^{\ast}A^{\ast}S^{\ast}S\right)$
$\displaystyle=AV_{m}\left(V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}\right)^{\dagger}\left(V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}\right)\left(V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}\right)^{\dagger}\left(V_{m}^{\ast}A^{\ast}S^{\ast}S\right).$

Due to the well-known property of the Moore–Penrose pseudoinverse that
$\left(V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}\right)\left(V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}\right)^{\dagger}$
is a projector onto
$\mathrm{Range}\left(V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}\right)$, it follows
that

$\displaystyle\left(V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}\right)^{\dagger}\left(V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}\right)\left(V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}\right)^{\dagger}=\left(V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}\right)^{\dagger}.$

Thus we have

$\displaystyle\widehat{\Phi}_{\mathcal{K}_{m}}^{2}$
$\displaystyle=AV_{m}\left(V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}\right)^{\dagger}\left(V_{m}^{\ast}A^{\ast}S^{\ast}S\right)=\widehat{\Phi}_{\mathcal{K}_{m}}.$

Similarly for $\widehat{\Pi}_{\mathcal{K}_{m}}$,

$\displaystyle\widehat{\Pi}_{\mathcal{K}_{m}}^{2}$
$\displaystyle=V_{m}\left(V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}\right)^{\dagger}\left(V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}\right)\left(V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}\right)^{\dagger}\left(V_{m}^{\ast}A^{\ast}S^{\ast}S\right)A$
$\displaystyle=V_{m}\left(V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}\right)^{\dagger}\left(V_{m}^{\ast}A^{\ast}S^{\ast}S\right)A=\widehat{\Pi}_{\mathcal{K}_{m}}.$

The range and null space containment assertions follow directly from the
definitions of $\widehat{\Pi}_{\mathcal{K}_{m}}$ and
$\widehat{\Phi}_{\mathcal{K}_{m}}$. However, we note that if
$\left(V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}\right)^{\dagger}$ is not a true
inverse, then it has as its null space
$\mathrm{Range}\left(V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}\right)^{\perp}$.
Furthermore, $S$ itself has a dimension $n-s$ null space. Thus, we can only
make null space containment statements, and these are not quite projectors.

We make one additional observation. Since
$\mathrm{Range}\left(AV_{m}\right)\subseteq\mathrm{Range}\left(B_{m+1}\right)$,
the sketching assumption holds for $\mathrm{Range}\left(AV_{m}\right)$, as
long as $\varepsilon$ is not too large. Thus, we prove the following.

######  Lemma 4.3.

For a sufficiently small sketching parameter $\varepsilon$ defined in (4), it
holds that $V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}$ is invertible.

######  Proof 4.4.

We use the notation $\mathbf{e}_{i}$ to denote the $i$-th Cartesian basis
vector and $AV_{m}\mathbf{e}_{i}$ to be the $i$-th column of $AV_{m}$. The
$(k,\ell)$-th entry of $V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}$ is of the form
$\langle SAV_{m}\mathbf{e}_{\ell},SAV_{m}\mathbf{e}_{k}\rangle$. From (3), it
follows

$\displaystyle\langle
SAV_{m}\mathbf{e}_{\ell},SAV_{m}\mathbf{e}_{k}\rangle=\langle
AV_{m}\mathbf{e}_{\ell},AV_{m}\mathbf{e}_{k}\rangle+\delta_{k\ell}\left\|AV_{m}\mathbf{e}_{k}\right\|\left\|AV_{m}\mathbf{e}_{\ell}\right\|,$

where for all $k$ and $\ell$, $|\delta_{k\ell}|<\varepsilon$. This is a
perturbation of a nonsingular Gram matrix, and for $\varepsilon$ sufficiently
small, the perturbation is still invertible.

Thus, while these idempotent matrices are not projectors, they do act as
projectors _when restricted_ to subspaces satisfying the sketching assumption
(4).

###  4.3 Semi-inner product induced by $S^{\ast}S$

How shall we understand the almost-projectors
$\widehat{\Pi}_{\mathcal{K}_{m}}$ and $\widehat{\Phi}_{\mathcal{K}_{m}}$?
Consider the fact that if $S^{\ast}S$ were nonsingular, it would induce an
inner product. In this case, one can easily show that these would be true
projectors with $\widehat{\Pi}_{\mathcal{K}_{m}}$ being an orthogonal
projector onto $\mathrm{Range}(V_{m})$ with respect to the
$A^{\ast}S^{\ast}SA$-inner-product, and $\widehat{\Phi}_{\mathcal{K}_{m}}$
being an orthogonal projector onto $\mathrm{Range}(AV_{m})$ with respect to
the $S^{\ast}S$-inner-product. This would simply be the method called often in
the literature _weighted GMRES_[11].

However, $S^{\ast}S$ is _necessarily_ singular in the case of sketching. It is
thus pointed out in [3] that it induces a nonnegative bilinear form satisfying
all the axioms of inner products except that
$\mathbf{w}^{\ast}S^{\ast}S\mathbf{w}=0$ does not necessarily imply that
$\mathbf{w}=\mathbf{0}$. Of course, because of the nature of the sketching
process, this is not surprising. Thus, sketched GMRES may be understood as an
extreme version of weighted GMRES [11] wherein most of the weight values are
actually zero, and the nonzero weights have been chosen such that the weighted
semi-inner product still acts as an inner product on the correction and
constraint spaces.

It is observed in [3] that the bilinear form induced by positive semi-definite
$S^{\ast}S$ is called a _semi-inner product_ or pseudo-inner product. We could
then call the matrices $\widehat{\Pi}_{\mathcal{K}_{m}}$ and
$\widehat{\Phi}_{\mathcal{K}_{m}}$ _semi-projectors_ induced by semi-inner
products. In the case that $S$ is a sketching operator, we might also consider
the terminology _sketched projectors_. Let us denote the semi-inner product by
$\langle\mathbf{u},\mathbf{w}\rangle_{S^{\ast}S}:=\mathbf{w}^{\ast}S^{\ast}S\mathbf{u}$.
In [18], the authors point out further that the sketching assumption (4) means
that this semi-inner product acts as an inner product on the subspace in
question. We can interpret Lemma 4.3 as a consequence of this fact; a Gram
matrix formed from any set of basis vectors for which (4) holds, using the
sketched semi-inner product, is nonsingular.

###  4.4 Relationship between residual angle and sketched constraint space

In [10], the authors explored many aspects of the geometric properties of
quantities and subspaces generated by Krylov subspace iterations, in the
context of understanding aspects of convergence. We now extend some of these
results to the sketched GMRES setting. Much work is done concerning geometric
aspects of minimum residual methods such as GMRES, and we focus mainly on
extending the ideas used to develop [10, Theorem 4.5] and results subsequent
to that.

We denote by $\theta\left(\mathbf{v},\,\mathbf{u}\right)$ the angle between
vectors $\mathbf{v}$ and $\mathbf{u}$, and by
$\theta\left(\mathbf{v},\,\mathcal{W}\right)$ the principal angle between the
vector $\mathbf{v}$ and the subspace $\mathcal{W}$. The authors of [10] make
use of the observation (true for any orthogonal projector onto a subspace)
that

$\displaystyle\left\|\Phi_{\mathcal{K}_{m}}\mathbf{r}_{0}\right\|=\left\|\mathbf{r}_{0}\right\|\cos\theta\left(\mathbf{r}_{0},\,A\mathcal{K}_{m}\right)\
\mbox{and,}\
\left\|\left(I-\Phi_{\mathcal{K}_{m}}\right)\mathbf{r}_{0}\right\|=\left\|\mathbf{r}_{0}\right\|\sin\theta\left(\mathbf{r}_{0},\,A\mathcal{K}_{m}\right),$

This identity does not necessarily hold for sketched semi-projectors, but the
residual projector $\Phi_{\mathcal{K}_{m}}$ can be related to a true projector
acting in the sketched constraint space. This can be used to connect us back
to the theory used in [10, Section 4], when acting on vectors from the
sketched subspace, specifically $\mathbf{r}_{0}$.

######  Lemma 4.5.

For sketched GMRES, it holds for the sketched residual projection

(9)
$\displaystyle\left\|\widetilde{\mathbf{r}}_{m}\right\|_{S^{\ast}S}=\left\|\left(I-\widetilde{\Phi}_{\mathcal{K}_{m}}\right)\mathbf{r}_{0}\right\|_{S^{\ast}S}=\left\|S\mathbf{r}_{0}\right\|\sin\theta\left(S\mathbf{r}_{0},\,SA\mathcal{K}_{m}\right).$

######  Proof 4.6.

This follows directly from the fact that
$S\widetilde{\Phi}_{\mathcal{K}_{m}}=\Phi_{\widetilde{\mathcal{K}}_{m}}S$,
where $\Phi_{\widetilde{\mathcal{K}}_{m}}$ is the orthogonal projector acting
in the sketched space, projecting onto $SA\mathcal{K}_{m}(A,\mathbf{r}_{0})$.
From this the result follows directly since

$\displaystyle\left\|S\left(I-\widetilde{\Phi}_{\mathcal{K}_{m}}\right)\mathbf{r}_{0}\right\|=\left\|\left(I-\Phi_{\widetilde{\mathcal{K}}_{m}}\right)S\mathbf{r}_{0}\right\|,$

to which we can then apply (9).

Notice that we can consider developing bounds on the true residual in terms of
sketched angles. To make use of this, we first relate trigonometric functions
of sketched angles back to their unsketched counterparts.

######  Lemma 4.7.

Let $\mathbf{u}$ and $\mathbf{v}$ satisfy the sketching assumption (4) for
$S$. Then it holds that

$\displaystyle\dfrac{\cos\theta\left(\mathbf{u},\,\mathbf{v}\right)-\varepsilon}{\left(1+\varepsilon\right)^{2}}$
$\displaystyle\leq\cos\theta\left(S\mathbf{u},\,S\mathbf{v}\right)\leq\dfrac{\cos\theta\left(\mathbf{u},\,\mathbf{v}\right)+\varepsilon}{\left(1-\varepsilon\right)^{2}},\qquad\mbox{and}$
(10)
$\displaystyle\dfrac{\sqrt{\sin^{2}\theta\left(\mathbf{u},\,\mathbf{v}\right)-10\varepsilon}}{\left(1-\varepsilon\right)^{2}}$
$\displaystyle\leq\sin\theta\left(S\mathbf{u},\,S\mathbf{v}\right)\leq\dfrac{\sqrt{\sin^{2}\theta\left(\mathbf{u},\,\mathbf{v}\right)+16\varepsilon}}{\left(1+\varepsilon\right)^{2}}.$

######  Proof 4.8.

Dividing the inequality (3) by
$\left\|\mathbf{u}\right\|\left\|\mathbf{v}\right\|$ yields

(11)
$\displaystyle\cos\theta\left(\mathbf{u},\,\mathbf{v}\right)-\varepsilon\leq\dfrac{\left\|S\mathbf{u}\right\|\left\|S\mathbf{v}\right\|}{\left\|\mathbf{u}\right\|\left\|\mathbf{v}\right\|}\cos\theta\left(S\mathbf{u},\,S\mathbf{v}\right)\leq\cos\theta\left(\mathbf{u},\,\mathbf{v}\right)+\varepsilon.$

From the sketching assumption (4), we have the bound

(12)
$\displaystyle\dfrac{1}{\left(1+\varepsilon\right)^{2}}\leq\dfrac{\left\|\mathbf{u}\right\|\left\|\mathbf{v}\right\|}{\left\|S\mathbf{u}\right\|\left\|S\mathbf{v}\right\|}\leq\dfrac{1}{\left(1-\varepsilon\right)^{2}}$

Multiplying (11) through by
$\dfrac{\left\|\mathbf{u}\right\|\left\|\mathbf{v}\right\|}{\left\|S\mathbf{u}\right\|\left\|S\mathbf{v}\right\|}$
and applying (12) yields the first result in (10).

To obtain the second result, we take the first result, square it, negate it
(reversing inequality order) and add $1$ to it, yielding

$\displaystyle
1-\left(\dfrac{\cos\theta\left(\mathbf{u},\,\mathbf{v}\right)+\varepsilon}{\left(1-\varepsilon\right)^{2}}\right)^{2}\leq\sin^{2}\theta\left(S\mathbf{u},\,S\mathbf{v}\right)\leq
1-\left(\dfrac{\cos\theta\left(\mathbf{u},\,\mathbf{v}\right)-\varepsilon}{\left(1+\varepsilon\right)^{2}}\right)^{2}.$

Getting everything over a common denominator and multiplying the numerators
out yields

(13)
$\displaystyle\dfrac{\sin^{2}\theta\left(\mathbf{u},\,\mathbf{v}\right)-4\varepsilon+7\varepsilon^{2}-4\varepsilon^{3}+\varepsilon^{4}-2\varepsilon\cos\theta\left(\mathbf{u},\,\mathbf{v}\right)}{\left(1-\varepsilon\right)^{4}}\leq\sin^{2}\theta\left(S\mathbf{u},\,S\mathbf{v}\right)\cdots$
(14)
$\displaystyle\cdots\leq\dfrac{\sin^{2}\theta\left(\mathbf{u},\,\mathbf{v}\right)+4\varepsilon+5\varepsilon^{2}+4\varepsilon^{3}+\varepsilon^{4}+2\varepsilon\cos\theta\left(\mathbf{u},\,\mathbf{v}\right)}{\left(1+\varepsilon\right)^{4}}$

As we assume $0<\varepsilon\ll 1$, we can write
$4\varepsilon+5\varepsilon^{2}+4\varepsilon^{3}+\varepsilon^{4}+2\varepsilon\cos\theta\left(\mathbf{u},\,\mathbf{v}\right)<16\varepsilon$,
and
$-4\varepsilon+7\varepsilon^{2}-4\varepsilon^{3}+\varepsilon^{4}-2\varepsilon\cos\theta\left(\mathbf{u},\,\mathbf{v}\right)\geq-10\varepsilon$.
Inserting these inequalities into (14) and taking square roots yields the
second result from (10).

Lemma 4.7 directly leads to a bound on the sketched GMRES residual norm.

######  Corollary 4.9.

The sketched GMRES residual $\widetilde{\mathbf{r}}_{m}$ satisfies the
subspace-angle-based bound

$\displaystyle\left\|\widetilde{\mathbf{r}}_{m}\right\|\leq\left\|\mathbf{r}_{0}\right\|\dfrac{\sqrt{\sin\theta\left(\mathbf{r}_{0},\,A\mathcal{K}_{m}\left(A,\mathbf{r}_{0}\right)\right)+16\varepsilon}}{1-\varepsilon^{2}}.$

######  Proof 4.10.

Combining (9) with the left side of the sketching assumption (4) for
$\widetilde{\mathbf{r}}_{m}$ yields

$\displaystyle\left(1-\varepsilon\right)\left\|\widetilde{\mathbf{r}}_{m}\right\|\leq\left\|S\mathbf{r}_{0}\right\|\sin\theta\left(S\mathbf{r}_{0},\,SA\mathcal{K}_{m}\left(A,\mathbf{r}_{0}\right)\right).$

Dividing through by $1-\varepsilon$ and applying the right side of the second
inequality in (10) yields the result.

###  4.5 Analysis based on the sketched Arnoldi relation

We note that our basis $V_{m}$ for $\mathcal{K}_{m}(A,\mathbf{r}_{0})$
satisfies

(15) $\displaystyle SAV_{m}=SV_{m+1}\underline{H_{m}}$

where $\underline{H_{m}}$ is constructed progressively using a truncated
Arnoldi orthogonalization with truncation parameter $t$; thus, is $t$-banded.
Computing the economy QR decomposition $SV_{m+1}=\widehat{W}_{m+1}R_{m+1}$ and
setting $\underline{\widehat{H}_{m}}=R_{m+1}\underline{H_{m}}$ (which is full
upper Hessenberg) yields the sketched Arnoldi relation

(16) $\displaystyle SAV_{m}=\widehat{W}_{m+1}\underline{\widehat{H}_{m}}.$

This is what is advocated for in, e.g., [2] but is not the approach we take in
this manuscript. However, we note that sketched GMRES algorithms derived from
(16) and (15) are mathematically equivalent. Thus, we analyze the formulation
of sketched GMRES derived from (16).

Consider that the version of sketched GMRES derived from (16) minimizes the
same functional,

(17) $\displaystyle\widetilde{\mathbf{x}}_{m}$
$\displaystyle=\underset{\mathbf{t}\in\mathcal{K}_{m}(A,\mathbf{r}_{0})}{\mathrm{argmin}}\left\|S\left(\mathbf{r}_{0}-A\mathbf{t}\right)\right\|$
$\displaystyle\iff\widetilde{\mathbf{x}}_{m}=V_{m}\widetilde{\mathbf{y}}_{m}\
\ \mbox{where}\ \ \widetilde{\mathbf{y}}_{m}$
$\displaystyle=\underset{\mathbf{y}\in\mathbb{C}^{m}}{\mathrm{argmin}}\left\|\underline{\widehat{H}_{m}}\mathbf{y}-\left\|S\mathbf{r}_{0}\right\|\mathbf{e}_{1}\right\|.$

The functional that $\widetilde{\mathbf{y}}_{m}$ minimizes is the standard
GMRES functional, simply for the sketched problem. _This is important_ because
any analysis of GMRES with proofs proceeding from an analysis of the
minimization of this functional can be mapped onto analogous results for
sketched GMRES. As with full GMRES, the minimization in (4.5) can be solved by
progressively QR-factorizing $\underline{\widehat{H}_{m}}$ using one Givens
rotation per iteration, obtaining at step $m$, yielding
$\underline{\widehat{H}_{m}}={\underline{\widehat{Q}}_{m}}{\underline{\widehat{R}_{m}}}$.

######  Theorem 4.11.

Let us write
${\underline{\widehat{Q}_{m}}}=\left(\underline{\hat{q}}_{ij}\right)\in\mathbb{C}^{(m+1)\times(m+1)}$.
Then it follows that

$\displaystyle\sin\theta\left(S\mathbf{r}_{0},\,SA\mathcal{K}_{m}(A,\mathbf{r}_{0})\right)=\left|q_{m+1,1}\right|.$

Furthermore, it follows that

$\displaystyle\sin\theta\left(S\mathbf{r}_{m-1},\,SA\mathcal{K}_{m}(A,\mathbf{r}_{0})\right)=s_{m}$

where $s_{m}$ is the $m$-th Givens sine used for constructing the QR-
factorization of $\underline{\widehat{H}_{m}}$. The sketched residual norm
satisfies
$\left\|S\tilde{\mathbf{r}}_{m}\right\|=|s_{m}|\left\|S\tilde{\mathbf{r}}_{m-1}\right\|$.

######  Proof 4.12.

We note that

$\displaystyle\theta\left(S\mathbf{r}_{0},\,SA\mathcal{K}_{m}(A,\mathbf{r}_{0})\right)=\theta\left(\mathbf{e}_{1},\,\mathrm{Range}\left(\underline{\widehat{H}_{m}}\right)\right).$

The rest follows step-for-step the proof of [10, Theorem 4.5] and analysis
thereafter.

Thus, if we implement sketched GMRES using (16) following [2], we can use the
Givens sine $s_{m}$ at each iteration $m$ to estimate the residual using
bounds such as in Lemma 4.7 and Corollary 4.9.

######  Corollary 4.13.

For the implementation of sketched GMRES using (16) following [2], the true
residual norm admits a cheaply computable estimate satisfying the progressive
updating scheme using the Givens sines, namely

(18)
$\displaystyle\left\|\tilde{\mathbf{r}}_{m}\right\|\leq\dfrac{s_{m}+4\varepsilon}{1-\varepsilon^{2}}\left\|\tilde{\mathbf{r}}_{m-1}\right\|.$

######  Proof 4.14.

This estimate follows from combining the results in Theorem 4.11 with the
sketching assumption (4) and the inequality relating sines of sketched angles
to sines of the unsketched analogs (10). Altogether, this yields

$\displaystyle\left\|\tilde{\mathbf{r}}_{m}\right\|\leq\dfrac{\sqrt{s_{m}^{2}+16\varepsilon}}{1-\varepsilon^{2}}\left\|\tilde{\mathbf{r}}_{m-1}\right\|.$

We then note that

$\displaystyle\sqrt{s_{m}^{2}+16\varepsilon}=\sqrt{\left(s_{m}+4\varepsilon\right)^{2}-8s_{m}\sqrt{\varepsilon}}\leq
s_{m}+4\varepsilon,$

with the last inequality being a tight bound for sufficiently small
$\varepsilon$. This yields the result.

We note that to use this bound requires knowledge of $\varepsilon$, which is
generally not know in the context of sketching. However, if we have a
practical upper bound for $\varepsilon$, we can use it to obtain a computable
update for the upper bound on $\left\|\mathbf{r}_{m}\right\|$. For example, if
we know that $\varepsilon<\varepsilon_{bound}$, it follows that

$\displaystyle\left\|\tilde{\mathbf{r}}_{m}\right\|\leq\dfrac{s_{m}+4\varepsilon}{1-\varepsilon^{2}}\left\|\tilde{\mathbf{r}}_{m-1}\right\|\leq\dfrac{s_{m}+4\cdot\varepsilon_{bound}}{1-\varepsilon_{bound}^{2}}\left\|\tilde{\mathbf{r}}_{m-1}\right\|$

This sort of relationship is useful in the case that we implement sketched
(augmented) GMRES using (16) following [2]. However, that is not what we
advocate in this paper. Thus, we must transfer this estimate to the version
using truncation.

An impetus in this context for using truncation is to keep the amount of work
done per iteration (e.g., inner products) constant rather than growing with
each iteration. Thus, any effort to exploit (18) should not increase the per-
iteration work in a way that grows with each iteration. Consider that to
obtain $\underline{\widehat{H}_{m}}$ would require that we run the full
sketched Arnoldi process, defeating the computational advantage we achieve by
truncating Arnoldi. Instead, we propose to use this estimate sparingly, as a
relatively cheap residual estimator to use when convergence is suspected based
on looser estimates. Consider that from the development of (15), we can obtain
$\underline{\widehat{H}_{m}}$ at any iteration by computing the QR-
factorization $SV_{m+1}=\widehat{W}_{m+1}R_{m+1}$ and obtaining
$\underline{\widehat{H}_{m}}=R_{m+1}\underline{H_{m}}$. The QR-factorization
of $\underline{\widehat{H}_{m}}$ can be computed by applying Givens rotations
to obtain the sequence of Givens sines, thereby enabling us to exploit
Corollary 4.13 to estimate the residual (using an estimate of $\varepsilon$).
If it is determined that the iteration has not converged,
$\underline{\widehat{H}_{m}}$, $\widehat{W}_{m+1}$, the Givens rotations and
can be stored for use in residual estimation for future iterations, to avoid
redundant work.

###  4.6 The relationship between sketched GMRES and sketched FOM

Interestingly, it also follows from (16) to derive sketched GMRES that we can
show that the well-known relationship between true GMRES and true FOM can be
extended to the sketched counterparts. However, the result is actually true
for any mathematically equivalent implementations of sketched GMRES and
sketched FOM [7].

For clarity, we identify quantities associated with sketched GMRES using the
superscript $\cdot^{(G)}$ and quantities associated with sketched FOM using
the superscript $\cdot^{(F)}$. By sketched FOM, we simply mean that we apply a
sketched version of the FOM constraint to define the FOM iterate
$\tilde{\mathbf{x}}_{m}^{(F)}=\mathbf{x}_{0}+\tilde{\mathbf{t}}_{m}^{(F)}$ via

$\displaystyle\mbox{Select}\ \
\tilde{\mathbf{t}}_{m}^{(F)}\in\mathcal{K}(A,\mathbf{r}_{0})\ \ \mbox{such
that}\ \
S\left(\mathbf{b}-A\left(\mathbf{x}_{0}+\tilde{\mathbf{t}}_{m}^{(F)}\right)\right)\perp
S\mathcal{K}_{m}(A,\mathbf{r}_{0}).$

For $\tilde{\mathbf{t}}_{m}^{(F)}=V_{m}\tilde{\mathbf{y}}_{m}^{(F)}$ where
$\tilde{\mathbf{y}}_{m}^{(F)}$ satisfies the FOM condition

$\displaystyle\widehat{H}_{m}\tilde{\mathbf{y}}_{m}^{(F)}=\left\|S\mathbf{r}_{0}\right\|\mathbf{e}_{1},$

where $\widehat{H}_{m}$ is the square matrix formed by the first $m$ rows of
${\underline{\widehat{H}_{m}}}$. This means that sketched GMRES and sketched
FOM form an MR/OR pair, as described in, e.g., [10, Section 1]. We omit
extending every such result, but we point the reader to the results in [6],
which are proven by comparing the QR-factorizations of the square and
rectangular upper Hessenberg matrices. In particular, following the derivation
of [22, eq. 6.74], a comparison of the QR-factorizations of
${\underline{\widehat{H}_{m}}}$ and $\widehat{H}_{m}$ shows the following.

######  Theorem 4.15.

The the sketched GMRES and sketched FOM iterates have the same relationship as
the true GMRES and FOM iterates, namely

$\displaystyle\tilde{\mathbf{x}}_{m}^{(G)}=c_{m}^{2}\tilde{\mathbf{x}}_{m}^{(F)}+s_{m}^{2}\tilde{\mathbf{x}}_{m-1}^{(G)}.$

##  5 Analysis of sketched augmented GMRES

We have established that the sketched GMRES method is not a projection method,
but it is a semi-projection method. In the context of the characterization of
augmented residual minimization methods from [25] that we summarize in Section
2.2, we may ask: does a sketched augmented minimization method fit into a
generalization of the framework from [25]?

###  5.1 Sketched augmented GMRES via projection framework

Consider the sketched augmented GMRES approximation
$\widetilde{\mathbf{x}}_{m}=\mathbf{x}_{0}+U\widetilde{\mathbf{z}}_{m}+V_{m}\widetilde{\mathbf{y}}_{m}$.
The vectors $\widetilde{\mathbf{z}}_{m}$ and $\widetilde{\mathbf{y}}_{m}$ are
components of a solution to the normal equations

(19)
$\displaystyle\begin{bmatrix}U^{\ast}A^{\ast}S^{\ast}SAU&U^{\ast}A^{\ast}S^{\ast}SAV_{m}\\\
V_{m}^{\ast}A^{\ast}S^{\ast}SAU&V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}\end{bmatrix}\begin{bmatrix}\widetilde{\mathbf{z}}\\\
\widetilde{\mathbf{y}}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}U^{\ast}A^{\ast}S^{\ast}S\mathbf{r}_{0}\\\
V_{m}^{\ast}A^{\ast}S^{\ast}S\mathbf{r}_{0}\end{bmatrix}.$

Because this is a sketched minimum residual method, these normal equations’
coefficient matrix may have a null space; thus there could be multiple
solutions. By choosing

$\displaystyle\begin{bmatrix}\widetilde{\mathbf{z}}_{m}\\\
\widetilde{\mathbf{y}}_{m}\end{bmatrix}=\begin{bmatrix}U^{\ast}A^{\ast}S^{\ast}SAU&U^{\ast}A^{\ast}S^{\ast}SAV_{m}\\\
V_{m}^{\ast}A^{\ast}S^{\ast}SAU&V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}\end{bmatrix}^{\dagger}\begin{bmatrix}U^{\ast}A^{\ast}S^{\ast}S\mathbf{r}_{0}\\\
V_{m}^{\ast}A^{\ast}S^{\ast}S\mathbf{r}_{0}\end{bmatrix},$

the Moore-Penrose pseudo-inverse solution, we obtain the solution which has
minimum norm, containing no null space components. However, this consideration
is moot when (4) holds.

######  Corollary 5.1 (to Lemma 4.3).

For sketching parameter $\varepsilon$ sufficiently small, if the sketching
assumption in (4) holds for $\mathrm{Range}\left(AU\right)$, then the
coefficient matrix in (19) is invertible, and so is the $(1,1)$ block
$U^{\ast}A^{\ast}S^{\ast}SAU$.

In [25], the authors follow derivation first pointed out in [21] which uses
block Gaussian elimination to eliminate $\widetilde{\mathbf{z}}_{m}$ from the
second block row of equations, yielding the result. In the case that the
matrix has a null space, general theory for when block Gaussian elimination
can be applied to develop a generalized Schur complement has been derived;
see, e.g., [8]. The case discussed in the present work is simpler.

######  Corollary 5.2.

For sketching parameter $\varepsilon$ sufficiently small, the coefficient
matrix in (19) admits a Schur complement via elimination of the $(1,2)$ block
using the block row operation

$\displaystyle\begin{bmatrix}I&\\\
-V_{m}^{\ast}A^{\ast}S^{\ast}SAU\left(U^{\ast}A^{\ast}S^{\ast}SAU\right)^{-1}&I\end{bmatrix}\begin{bmatrix}U^{\ast}A^{\ast}S^{\ast}SAU&U^{\ast}A^{\ast}S^{\ast}SAV_{m}\\\
V_{m}^{\ast}A^{\ast}S^{\ast}SAU&V_{m}^{\ast}A^{\ast}S^{\ast}SAV_{m}\end{bmatrix},$

leading to the equivalent system

(20)
$\displaystyle\begin{bmatrix}U^{\ast}A^{\ast}S^{\ast}SAU&U^{\ast}A^{\ast}S^{\ast}SAV_{m}\\\
\mathbf{0}&V_{m}^{\ast}A^{\ast}S^{\ast}S\left(I-\widetilde{\Phi}_{U}\right)AV_{m}.\end{bmatrix}\begin{bmatrix}\widetilde{\mathbf{z}}\\\
\widetilde{\mathbf{y}}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}U^{\ast}A^{\ast}S^{\ast}S\mathbf{r}_{0}\\\
V_{m}^{\ast}A^{\ast}S^{\ast}S\left(I-\widetilde{\Phi}_{U}\right)\mathbf{r}_{0}\end{bmatrix}.$

Thus, as in [25], we have isolated in the second block of equations the
variables $\widetilde{\mathbf{y}}$ associated to the space
$\mathrm{Range}\left(V_{m}\right)$ that grows at each iterations. We can prove
the following theorem.

######  Theorem 5.3.

For $U$ such that $\mathrm{Range}\left(AU\right)$ satisfies the sketching
assumption (4), and for $m$ not too large and $\varepsilon$ sufficiently
small, applying augmented sketched GMRES for the augmented space
$\mathrm{Range}\left(\mathbf{V}_{m}\right)$ is equivalent to applying a
sketched GMRES iteration to the equation

(21)
$\displaystyle\left(I-\widetilde{\Phi}_{\mathbf{U}}\right)A\mathbf{t}=\left(I-\widetilde{\Phi}_{\mathbf{U}}\right)\mathbf{r}_{0},$

selecting
$\mathbf{t}_{m}=V_{m}\widetilde{\mathbf{y}}_{m}\in\mathrm{Range}\left(V_{m}\right)$
according to

(22)
$\displaystyle\widetilde{\mathbf{y}}_{m}=\underset{\mathbf{y}\in\mathbb{C}^{m}}{\mathrm{argmin}}\left\|S\left(I-\widetilde{\Phi}_{\mathbf{U}}\right)\left(\mathbf{r}_{0}-AV_{m}\mathbf{y}\right)\right\|$

and constructing
$\widetilde{\mathbf{x}}_{m}=\mathbf{x}_{0}+\widetilde{\Pi}_{U}\boldsymbol{\eta}_{0}+(I-\widetilde{\Pi}_{U})\widetilde{\mathbf{t}}_{m}$.

######  Proof 5.4.

This has mostly already been proven in the above discussion. We observe only
that the second block of equations in (20) defines $\widetilde{\mathbf{y}}$ to
be the exactly the minimizer described in (22). Inserting
$\widetilde{\mathbf{y}}_{m}$ into the first block of equations yields

$\displaystyle\widetilde{\mathbf{z}}$
$\displaystyle=U^{\ast}A^{\ast}S^{\ast}S\mathbf{r}_{0}-U^{\ast}A^{\ast}S^{\ast}SAV_{m}\widetilde{\mathbf{y}}_{m}$
$\displaystyle\iff U\widetilde{\mathbf{z}}$
$\displaystyle=\widetilde{\Pi}_{U}\boldsymbol{\eta}_{0}-\widetilde{\Pi}_{U}\mathbf{t}_{m},$

completing the proof.

It follows directly that the full augmented sketched GMRES residual is the
same as that of the projected subproblem, meaning that the residual behavior
of this method is completely dictated by the properties of the subproblem and
the the space $\mathrm{Range}\left(V_{m}\right)$.

######  Corollary 5.5.

Let $U$ and $AU$ satisfy the same conditions as in Theorem 5.3. Then it
follows the the full sketched augmented GMRES residual
$\widetilde{\mathbf{r}}_{m}=\mathbf{b}-A\widetilde{\mathbf{x}}_{m}$ is the
same as the residual from approximately solving (21) via (22); i.e.,

$\displaystyle\widetilde{\mathbf{r}}_{m}=\left(I-\widetilde{\Phi}_{\mathbf{U}}\right)\left(\mathbf{r}_{0}-A\mathbf{t}_{m}\right)$

######  Proof 5.6.

We compute

$\displaystyle\widetilde{\mathbf{r}}_{m}$
$\displaystyle=\mathbf{b}-A\widetilde{\mathbf{x}}_{m}$
$\displaystyle=\mathbf{b}-A(\mathbf{x}_{0}+\widetilde{\Pi}_{U}\boldsymbol{\eta}_{0}+(I-\widetilde{\Pi}_{U})\widetilde{\mathbf{t}}_{m})$
$\displaystyle=\mathbf{r}_{0}-\widetilde{\Phi}_{U}\mathbf{r}_{0}-(I-\widetilde{\Phi}_{U})A\widetilde{\mathbf{t}}_{m}),$

where we have taken advantage of the fact that
$A\widetilde{\Pi}_{U}=\widetilde{\Phi}_{U}A$.

We make the observation that, were we to build $V_{m}$ such that

$\displaystyle\mathrm{Range}\left(V_{m}\right)=\mathcal{K}_{m}\left(\left(I-\widetilde{\Phi}_{\mathbf{U}}\right)A,\left(I-\widetilde{\Phi}_{\mathbf{U}}\right)\mathbf{r}_{0}\right),$

we could use Corollary 5.5 to exploit the results from Theorem 4.11 and
thereafter directly, as the method could directly be interpreted as a fully
sketched version of GCRO-DR. Corollary 5.5 then implies that the residual
convergence behavior of the iterative method is equivalent to that of sketched
GMRES applied to (3), and we can estimate the resdidual norm using the Givens
sines, as described in Theorem 4.11. _However, we do not advocate building a
sketched projected Krylov subspace._

We have demonstrated that we can gain computational advantage by forgoing the
building of a sketched projected Krylov subspace. Instead, we have

$\displaystyle\mathrm{Range}\left(V_{m}\right)=\mathcal{K}_{m}\left(A,\mathbf{r}_{0}\right).$

The price for this advantage is that Corollary 5.5 no longer implies that the
iterative method is equivalent to a sketched GMRES iteration applied to (3).
However, all is not lost; we can still use Theorem 4.11. Consider that via
truncated sketched Arnoldi, we have the relation

$\displaystyle
SA\mathbf{V}_{m}=\begin{bmatrix}SAU&SV_{m+1}\end{bmatrix}\begin{bmatrix}I&\\\
&\underline{H_{m}}\,\end{bmatrix}.$

We compute the economy QR-factorization
$\begin{bmatrix}SAU&SV_{m+1}\end{bmatrix}=\widehat{\mathbf{W}}_{m+1}R_{m}$,
and block the upper triangular matrix

$\displaystyle R_{m}=\begin{bmatrix}R_{m}^{(1,1)}&R_{m}^{(1,2)}\\\
&R_{m}^{(2,2)}\end{bmatrix},\qquad\mbox{with}\qquad
R_{m}^{(1,1)}\in\mathbb{C}^{k\times k}.$

Then we can construct the upper-Hessenberg matrix

$\displaystyle\underline{G_{m}}=R_{m}\begin{bmatrix}I&\\\
&\underline{H_{m}}\,\end{bmatrix}=\begin{bmatrix}R_{m}^{(1,1)}&R_{m}^{(1,2)}\underline{H_{m}}\\\
&R_{m}^{(2,2)}\underline{H_{m}}\,\end{bmatrix}=\begin{bmatrix}R_{m}^{(1,1)}&\ast\\\
&\underline{\widehat{H}_{m}}\end{bmatrix},$

where $\underline{\widehat{H}_{m}}:=R_{m}^{(2,2)}\underline{H_{m}}$. It
follows that we have the relation

$\displaystyle SA\mathbf{V}_{m}=\widehat{\mathbf{W}}_{m+1}\underline{G_{m}}.$

We observe that $\underline{G_{m}}$ is upper Hessenberg but with an already-
triangular (1,1)-block. This means that a Givens rotation-based QR
factorization of $\underline{G_{m}}$ reduces to applying Givens rotations to
$\underline{\widehat{H}_{m}}$. The proof of Theorem 4.11 allows us to conclude
that

$\displaystyle\theta\left(S\mathbf{r}_{0},\,\mathrm{Range}\left(SA\mathbf{V}_{m}\right)\right)=\theta\left(\mathbf{e}_{1},\,\mathrm{Range}\left(\underline{G_{m}}\right)\right)=\theta\left(\mathbf{e}_{1},\,\mathrm{Range}\left(\underline{Q_{m}}\right)\right),$

where $\underline{Q_{m}}\in\mathbb{C}^{(m+1)\times(m+1)}$ is the Q-factor of
the full QR factorization of $\underline{G_{m}}$. Thus, because of the
partially-triangularized structure of $\underline{G_{m}}$, we can use the same
Givens rotation estimator strategy shown in Corollary 4.13. We demonstrate the
effectiveness of the residual estimator for an example problem in Figure 2.

######  Remark 5.7.

Note that in Figure 2, we compute the estimate for every iteration for
demonstrative purposes. In reality, the frequency of the use of this estimate
would depend on which implementation of augmented sketched GMRES. If an
orthonormal basis for $\mathrm{Range}\left(SA\mathbf{V}_{m}\right)$ is
unavailable, the estimate would be used sparingly, due to the aforementioned
expense of obtaining it. However, we note that in Algorithm 1, a basis-
whitening technique is used which does produce an orthonormal basis. In that
case, Givens rotations could be computed progressively, one-per-iteration, in
the usual way for GMRES. This would allow the estimate to be computed
inexpensively.

It is also of interest that given some true residual norm computations, one
could use the bound from Corollary 4.13 to get an estimate for $\varepsilon$.
This is a possible alternative method to the method for estimation of
$\varepsilon$ advocated in Section 3.4, which is what is actually implemented
in our proposed algorithm.

Figure 2: Testing the Givens-rotation-based residual estimator
$\left\|\tilde{\mathbf{r}}_{m}\right\|\leq\frac{s_{m}+4\varepsilon}{1-\varepsilon^{2}}\left\|\tilde{\mathbf{r}}_{m-1}\right\|$
for a problem of Neumann type from the experiments in Section 6.2 for various
values of $\varepsilon$.

##  6 Numerical experiments

We now present numerical experiments demonstrating the effectiveness of GMRES-
SDR. We compare to MATLAB’s gmres, the RANDGMRES method in the randKrylov
package111https://github.com/obalabanov/randKrylov, version as of November
2023, GCRO-DR and GMRES-DR. All experiments were performed in MATLAB Online,
using version R2023B.

###  6.1 Stokes problem

In the first experiment we solve a single linear system
$A\mathbf{x}=\mathbf{b}$, where $A$ is the vas_stokes_1M matrix of size
$N=1,090,664$ from the SuiteSparse Matrix Collection [9]. We solve the system
with ILU preconditioning to a target relative residual of $10^{-6}$. The
residual curves are shown in Figure 3 (left). The true residual norms for
RANDGMRES and GMRES-SDR are only computed and plotted at the end of each
restart cycle as within each cycle only the sketched residual is cheaply
available.

Figure 3: Convergence curves obtained from solving a single linear system
$A\mathbf{x}=\mathbf{b}$ (left), and a sequence of $5$ linear systems (right)
with the vas_stokes_1M matrix, to a residual tolerance of $10^{-6}$. All non-
augmented methods use a maximum number of $m=100$ Arnoldi iterations, while
the augmented methods use an augmentation subspace of dimension $k=20$, and
take a maximum of $m-k$ iterations. All sketching methods take $s=10(m+k)$ and
an Arnoldi truncation parameter $t=2$. A maximum of $10$ restarts is allowed
for all methods.

It is clear from Figure 3 (left) that the convergence of GMRES-SDR is
accelerated by augmented restarting, and it produces convergence curves
comparable with GCRO-DR and GMRES-DR.

In Table 1 we record the total computational resources utilized by each method
to converge, which we measure using the total number of matrix-vector products
with $A$ (denoted MATVECS or MV), the total number of inner products (IP) of
size $N$, and the overall runtime (T) in seconds, required to produce the
convergence curves in Figure 3 (left).

| GCRO-DR | RANDGMRES | GCRO-DR | GMRES-DR | GMRES-SDR
---|---|---|---|---|---
MV | 1,000 | 1,000 | 670 | 674 | 656
IP | 55,550 | 11 | 51,138 | 45,229 | 1,955
T (s) | 287.6 | 408.8 | 196.6 | 190.7 | 178.9
Table 1: Performance metrics for solving a single linear system
$A\mathbf{x}=\mathbf{b}$ with the vas_stokes_1M matrix to a residual tolerance
of $10^{-6}$ and with ILU preconditioning.

We also solve a sequence of $10$ systems $A\mathbf{x}^{(i)}=\mathbf{b}^{(i)}$
with the same fixed matrix and randomly generated right-hand sides (unit
Gaussian). The performance metrics are given in Table 2.

| GMRES | RANDGMRES | GCRO-DR | GMRES-DR | GMRES-SDR
---|---|---|---|---|---
MV | 10,000 | 10,000 | 2,646 | 6,943 | 3,439
IP | 555,500 | 110 | 208,398 | 468,478 | 10,255
T (s) | 2535.2 | 4070.8 | 1474.7 | 1924.4 | 914.8
Table 2: Performance metrics for solving a sequence of $10$ linear systems
$A\mathbf{x}^{(i)}=\mathbf{b}^{(i)}$ with the vas_stokes_1M matrix to a
residual tolerance of $10^{-6}$ and with ILU preconditioning.

From Tables 1 and 2 we see in that the combination of recycling and sketching
in GMRES-SDR leads to a beneficial reduction in MATVECS and inner products,
allowing for the fastest runtime for a single system or a sequence of systems.

In Figure 3 (right) we plot the relative residual curves of both GCRO-DR and
GMRES-SDR for the first $5$ systems in the sequence of Stokes problems. We see
that the convergence of both methods is almost identical, and it is
significantly improved with recycling beyond the first problem in the
sequence.

###  6.2 Neumann problem

We now solve $50$ linear systems $A\mathbf{x}^{(i)}=\mathbf{b}^{(i)}$ for
fixed matrix $A$ constructed as $A=D+cI$ where $D$ is the “Neumann” matrix of
size $10,609$ taken from MATLAB’s gallery, $I$ is the identity matrix and
$c=0.0001$. The right-hand side vectors $\mathbf{b}^{(i)}$ were generated with
random unit Gaussian entries. The utilized computational resources are shown
in Table 3. The reported timings are averages over $10$ runs of each solver.

| GMRES | RANDGMRES | GCRO-DR | GMRES-DR | GMRES-SDR
---|---|---|---|---|---
MV | 50,000 | 50,000 | 8,978 | 20,528 | 6,906
IP | 2,777,500 | 550 | 531,179 | 1,164,290 | 20,556
T (s) | 122.1 | 74.7 | 17.3 | 34.2 | 11.8
Table 3: Performance metrics for solving $50$ linear systems
$A\mathbf{x}^{(i)}=\mathbf{b}^{(i)}$ with the Neumann matrix to a residual
tolerance of $10^{-6}$. All non-augmented methods use a maximum number of
$m=100$ Arnoldi iterations, while the augmented methods use an augmentation
subspace of dimension $k=20$, and take a maximum of $m-k=80$ iterations. All
sketching methods use $s=10(m+k)$ and an Arnoldi truncation parameter $t=2$. A
maximum of $10$ restarts is allowed for all methods.

In Figure 4 we plot the relative residual curves obtained from solving both
the first (left) and last (right) problem, respectively. We see that the
GMRES-SDR convergence behaves in the same way as GCRO-DR, and substantially
improves by the time the last problem is solved. We observe that GCRO-DR
converges in slightly fewer iterations than GMRES-SDR, although as shown in
Table 3, the overall runtime of GMRES-SDR is considerably lower. We also note
that while GMRES-DR improves the convergence through deflated restarts for
each problem separately, GCRO-DR and GMRES-SDR lead to an overall faster
convergence as they benefit from a recycling subspace that is continuously
improved as the problem sequence progresses. This is a key characteristic of a
successful recycling method.

Figure 4: Convergence curves for the first (left) and last (right) Neumann
problem.

###  6.3 Convection–diffusion problem

Finally, we consider a problem with varying system matrix $A^{(i)}$ arising
from the finite-difference discretization of a 2D convection–diffusion problem
for different convection strengths. More precisely, we use $A^{(i)}=(L\otimes
I+I\otimes L)+\alpha^{(i)}(D\otimes I+I\otimes D)$, where

$L=(n+1)^{2}\begin{bmatrix}-2&1\\\ 1&-2&\ddots\\\ &\ddots&\ddots&1\\\
&&1&-2\end{bmatrix},\ D=\frac{n+1}{2}\begin{bmatrix}0&1\\\ -1&0&\ddots\\\
&\ddots&\ddots&1\\\ &&-1&0\end{bmatrix}\in\mathbb{R}^{n\times n},\ n=500.$

The vector $\mathbf{b}$ is chosen as the vector of all ones. The convection
strength is chosen as $\alpha^{(0)}=0,\alpha^{(1)}=5,\alpha^{(2)}=20$. The
convergence of GCRO-DR and GMRES-SDR is plotted in Figure 5 and the
performance metrics are listed in Table 4. In all cases, we a use Krylov
dimension of $m=80$ and an augmentation dimension of $k=20$. For GMRES-SDR we
consider the two variants discussed in Section 3.6. The first one, labeled
“inexact”, corresponds to using $U,SU,SA^{(i)}U$ as the augmentation for
problem $i+1$. In other words, we are _not_ computing $SA^{(i+1)}U$ with the
new matrix $A^{(i+1)}$.

There are several interesting observations to be made. First of all, GCRO-DR
converges robustly for all three problems, and the convergence seems to indeed
benefit from the augmentation. For the first two problems with $A^{(1)}$ and
$A^{(2)}$, GMRES-SDR (inexact) and GMRES-SDR (exact) converge almost
identically and are comparable to GCRO-DR. However, for the third problem with
$A^{(3)}$, GMRES-SDR (inexact) fails with an erratically increasing residual.
Note how this increase is not tracked by the residuals of the sketched
problems (shown as the yellow dotted curve), which are still decreasing within
each restart cycle. As explained in Section 3.6, the inexactness destroys the
relation (7) between the true and sketched residual problems. This is not the
case for GMRES-SDR (exact), which decreases the sketched residual steadily
(and the true residual almost monotonically) as expected from a GMRES-type
method. GMRES-SDR (inexact) fails because the change in going from $A^{(2)}$
to $A^{(3)}$ is “too large”, while it appears to be okay to go from $A^{(1)}$
to $A^{(2)}$ without computing $SA^{(2)}U$ explicitly.

Figure 5: Convergence curves for three convection–diffusion problems with increasing convection strength. For the two GMRES-SDR variants, there are two curves for each problem: (i) a dashed curve with “*” markers showing the true residual after each restart cycle, and (ii) a dotted curve corresponding to the residual of the sketched problem. | GCRO-DR | GMRES-SDR | GMRES-SDR
---|---|---|---
|  | (inexact) | (exact)
MATVECS | 1,345 | 1,886 | 1,971
Inner products | 105,695 | 5,613 | 4,668
Time (s) | 38.7 | 30.8 | 27.0
Table 4: Solving three convection–diffusion problems
$A^{(i)}\mathbf{x}^{(i)}=\mathbf{b}$. GMRES-SDR (inexact) fails to converge to
the target residual norm of $10^{-2}$ for the third problem.

##  7 Conclusions and future work

We have introduced GMRES-SDR, a GMRES variant that combines sketching with
deflated restarting. Our numerical tests indicate that GMRES-SDR can improve
over GCRO-DR and GMRES-DR in terms of arithmetic cost and runtime. The runtime
reduction compared to GCRO-DR (the best competing method) ranged between 30%
and 38% on the problems we considered. This is primarily the result of a
reduction in inner products due to sketching, and sometimes even a reduction
in matrix-vector products due to improved deflation (see, e.g., Table 3 where
GMRES-SDR required about 23% fewer matrix-vector products than GCRO-DR).

Our analysis extended and generalized several results known for GMRES to the
sketched and augmented case. In particular, we characterized GMRES-SDR as a
projection method using a semi-inner product.

In future work we would like to understand better the behaviour of the
“inexact” GMRES-SDR variant. In particular, it would be useful to know a
priori what constitutes a sufficiently slow variation in $A^{(i)}$ so that the
inexact variant still converges.

## References

  * [1] O. Balabanov and L. Grigori.  Randomized block Gram-Schmidt process for solution of linear systems and eigenvalue problems.  Technical Report arXiv:2111.14641, 2021.
  * [2] O. Balabanov and L. Grigori.  Randomized Gram-Schmidt process with application to GMRES.  SIAM J. Sci. Comput., 44(3):A1450–A1474, 2022.
  * [3] O. Balabanov and A. Nouy.  Randomized linear algebra for model reduction. Part I: Galerkin methods and error estimation.  Adv. Comput. Math., 45:2969–3019, 2019.
  * [4] O. Balabanov and A. Nouy.  Randomized linear algebra for model reduction–part II: minimal residual methods and dictionary-based approximation.  Adv. Comput. Math., 47:1–54, 2021.
  * [5] M. Bolten, E. de Sturler, C. Hahn, and M. L. Parks.  Krylov subspace recycling for evolving structures.  Comput. Methods Appl. Mech. Engrg., 391:114222, 2022.
  * [6] P. Brown.  A theoretical comparison of the Arnoldi and GMRES algorithms.  SIAM J. Sci. Stat. Comput., 12(1):58–78, 1991.
  * [7] L. Burke and S. Güttel.  Krylov subspace recycling with randomized sketching for matrix functions.  arXiv EPrint arXiv:2308.02290, arXiv, 2023.
  * [8] D. Carlson, E. Haynsworth, and T. Markham.  A generalization of the Schur complement by means of the Moore–Penrose inverse.  SIAM J. Appl. Math., 26(1):169–175, 1974.
  * [9] T. A. Davis and Y. Hu.  The University of Florida sparse matrix collection.  ACM Trans. Math. Software, 38(1):1–25, 2011.
  * [10] M. Eiermann and O. G. Ernst.  Geometric aspects of the theory of Krylov subspace methods.  Acta Numer., 10:251–312, 2001.
  * [11] A. Essai.  Weighted FOM and GMRES for solving nonsymmetric linear systems.  Numerical Algorithms, 18(3-4):277–292, 1998.
  * [12] A. Gaul.  Recycling Krylov subspace methods for sequences of linear systems.  PhD thesis, Technische Universität Berlin, Fakultät II - Mathematik und Naturwissenschaften, 2014.
  * [13] S. Güttel and M. Schweitzer.  Randomized sketching for Krylov approximations of large-scale matrix functions.  SIAM J. Matrix Anal. Appl., 44(3):1073–1095, 2023.
  * [14] P.-G. Martinsson and J. A. Tropp.  Randomized numerical linear algebra: Foundations and algorithms.  Acta Numer., 29:403–572, 2020.
  * [15] R. B. Morgan.  GMRES with deflated restarting.  SIAM J. Sci. Comput., 24(1):20–37, 2002.
  * [16] Y. Nakatsukasa and J. A. Tropp.  Fast & accurate randomized algorithms for linear systems and eigenvalue problems.  Technical Report arXiv:2111.00113, 2022.
  * [17] C. C. Paige, B. N. Parlett, and H. A. van der Vorst.  Approximate solutions and eigenvalue bounds from Krylov subspaces.  Numer. Linear Algebra Appl., 2(2):115–133, 1995.
  * [18] D. Palitta, M. Schweitzer, and V. Simoncini.  Sketched and truncated polynomial Krylov methods: Evaluation of matrix functions.  Technical Report arXiv:2306.06481, 2023.
  * [19] M. L. Parks.  The iterative solution of a sequence of linear systems arising from nonlinear finite element analysis.  University of Illinois at Urbana-Champaign, 2005.
  * [20] M. L. Parks, E. de Sturler, G. Mackey, D. D. Johnson, and S. Maiti.  Recycling Krylov subspaces for sequences of linear systems.  SIAM J. Sci. Comput., 28(5):1651–1674, 2006.
  * [21] M. L. Parks, K. M. Soodhalter, and D. B. Szyld.  A block recycled GMRES method with investigations into aspects of solver performance.  arXiv preprint arXiv:1604.01713, 2016.
  * [22] Y. Saad.  Iterative Methods for Sparse Linear Systems, 2nd edition.  SIAM, Philadelphia, 2000.
  * [23] Y. Saad and M. Schultz.  GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems.  SIAM J. Sci. Stat. Comput., 7(3):856–869, 1986.
  * [24] T. 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.
  * [25] K. M. Soodhalter, E. de Sturler, and M. E. Kilmer.  A survey of subspace recycling iterative methods.  GAMM-Mitt., 43(4):e202000016, 29, 2020.
  * [26] D. P. Woodruff.  Sketching as a tool for numerical linear algebra.  Found. Trends Theor. Comput. Sci., 10(1–2):1–157, 2014.


# WEKA-Based: Key Features and Classifier for French of Five Countries

Zeqian Li, Keyu Qiu, Chenxu Jiao, Wen Zhu, Haoran Tang

###### Abstract

This paper describes a French dialect recognition system that will
appropriately distinguish between different regional French dialects. A corpus
of five regions - Monaco, French-speaking, Belgium, French-speaking
Switzerland, French-speaking Canada and France, which is targeted for
construction by the Sketch Engine. The content of the corpus is related to the
four themes of eating, drinking, sleeping and living, which are closely linked
to popular life. The experimental results were obtained through the processing
of a python coded pre-processor and Waikato Environment for Knowledge Analysis
(WEKA) data analytic tool which contains many filters and classifiers for
machine learning.

##  I Introduction

With the development and popularity of artificial intelligence, dialect
recognition has become an important text task in the field of data mining[1].
There are many French dialects spoken around the World; French dialectologists
have studied varieties of local variations, but generally agree these cluster
into five main regional dialect: Monaco, Belgium, French-speaking areas in
both Switzerland and Canada and France. Since the spatial divide between these
countries, the development of French-speaking cultures in these regions has
been markedly different. These dialectal differences are worth discovering and
studying, not only in terms of cultural knowledge, but also in terms of their
importance in the field of computing. For example, dialect recognition can
help artificial intelligence to communicate effectively with users from
different geographical regions. To find these dialects’ unique features, we
use sketch engine to construct dialect corpus and use WEKA’s
StringToWordVector filter to extract words. What’s more, we used seven classic
classifiers — NaiveBayes, SMO, J48, ZeroR, JRip, c, and NaiveBayesMultinomial
to build the model. The use of WEKA’s powerful features has led to a series of
precise and valuable experimental results. Finally, advantages and performance
differences of each classifier are summarized through scientific and
quantitative assessment criteria for outcome data.

##  II Business Understanding

The objective of this task is to find the characteristics of different
regional French dialects to help people identify them accurately. The
requirement is to build a corpus of French dialects that meets the
requirements by means of the sketch engine. After obtaining the text
corresponding to the corpus, to convert it into aircraft rescue and
firefighting data file (.arff). Then, using the relevant methods and functions
of the WEKA data analytic tool to analyses classified French dialects corpus.
Additionally, comparing the performance of different classifiers through the
results. The definition of this data mining problem is How to effectively
identify the different dialects of French in different regions and help people
to be able to accurately determine their differences

##  III Data Preparation

Since Sketch Engine could only save .txt files, which can not be used directly
in WEKA. We use Python to clean these data, and process them into .arff
format.

The processing procedures and standard are as follows:

  * •

Read the .txt file by lines using Python, and filter out text without
or

.

  * •

Delete the
and
labels for each sentence.

  * •

Assume that a completed sentence is ended with a full point, exclamation mark
or question mark. According to this rule, filter all the unmatched sentence
again.

  * •

Delete all the sentences with repetition or too few words.

  * •

Use regex to match french characters, remove other symbols.

  * •

Add corresponding country label and quotation marks for each sentence.

  * •

Add header, composed of attributes, classification, and relationships of data
into.arff file applying the format of WEKA data set.

Thus, a dataset with 10052 instances was constructed. The number of instances
was 1674, 1011, 1588, 1478, 1143 for labels Morocco, Belgium, France, Canada,
and Switzerland respectively. Moreover, among these data, the maximum words in
a sentence is 1080, while the mean value is 40, and StdDev is 59. As for the
test set, The number of instances was 352, 110, 318, 284, 286 for labels with
same order as previous. The data statics of training data is shown in the
Table 1.

No. |  Label |  Count
---|---|---
1 |  MC |  1674
2 |  BE |  1011
3 |  FR |  1588
4 |  CA |  1478
5 |  CH |  1143
TABLE I: Training Data Statistic

##  IV Method

###  IV-A Filter

WEKA’s StringToWordVector filter was used in this experiment. It has the
ability to extract words from the document, and each word was given an
attribute of occurrence times.

###  IV-B Modeling

Based on WEKA’s classifier function and Alrehaili and Atwell’s (2017) research
on Arabic dialects, seven suitable classifiers, NaiveBayes, SMO, J48, ZeroR,
JRip, c, and NaiveBayesMultinomial, were selected for modeling, which
effectively avoided blind modeling and thus reduced the experimental time, the
result is shown in Table 2.

  * •

NaiveBayes: Naive Bayes method is a classification method based on Bayes
theorem and the independent assumption of feature conditions, which assumes
that each feature of samples is unrelated to other features.

  * •

SMO: SMO classifier implements the sequential minimum optimization algorithm,
which can be used to train support vector classifiers. This implementation
replaces all missing values globally and converts nominal attributes to binary
attributes [2].

  * •

J48: J48 classifier implements the C4.5 algorithm for building decision trees
which were pruned or unpruned. A pruned C4.5 decision tree was chosen in this
experiment.

  * •

ZeroR: ZeroR classifier relies on the target and ignores all predictions. It
only predicts most classes and has no predictability.

  * •

JRip JRip classifier implements a propositional rule learner, Repeated
Incremental Pruning to Produce Error Reduction (RIPPER) [3].

  * •

REPTree: REPTree classifier is a fast decision tree learner. It trims the
decision tree using error reduction pruning.

  * •

NaiveBayesMultinomial: NaiveBayesMultinomial is a classifier for building and
using a multinomial Naive Bayes classifier.

Classifier |  |  Evaluate on
---
Training Set
|  10 Fold Cross
---
Validation
|  60% Train
---
40% Test
NaiveBayes |  37.0683 |  35.6849 |  33.9229
SMO |  89.8607 |  59.2882 |  57.0222
J48 |  76.8297 |  49.6990 |  45.0106
ZeroR |  25.1875 |  25.1875 |  25.0264
JRip |  37.2162 |  34.2592 |  33.8173
REPTree |  58.6651 |  48.1571 |  44.7994
NaiveBayesMultinomial |  65.1178 |  59.0981 |  57.8669
TABLE II: The accuracy of different classifiers

###  IV-C Evaluating the performance of models

In this experiment, three validation method provided by WEKA were used to
evaluate the performance of the eight classifiers: use the training dataset,
10-fold cross-validation, and percentage split which divides the training set
into 60% for training and 40% for testing.

  * •

Evaluating on training set: A method that the training set and the test set
use the same data.

  * •

Cross-Validation: Divide the training set into N parts, use n-1 part for
training, use 1 part for testing, and so on for N times, and finally calculate
the overall result. N was set to 10 in this experiment.

  * •

Percentage Split: Divide the training set into two parts in proportion, one
for training and one for testing. In this experiment, the training set was
divided into 60% and 40% where 60% for training and 40% for testing.

###  IV-D Feature

The first experiment which is to choose the best classifier to recognize
French dialects show that SMO is the best machine learning classification
algorithm, but we could still improve the accuracy by adjusting the parameters
and features. The WordTokenizer setting assumes that the feature is strings or
words between spaces while CharacterNGramTokenizer assumes features are 1, 2
or 3 character sequences. NGramTokenizer tokenize the inputs into different
kinds of n-grams. The input can be a character vector of any length. We tried
to set Max and Min both to 1 gives a model based on single characters; Max and
Min both to 2 gives us a bigram model; Max and Min both to 3 is a trigram
model. Moreover, we also tried to use WordTokenizer with TF and IDF both False
and NGramTokenizer with Max to 3 and Min to 1.

The result of different tokenizer is shown in the Table 3.

Features |  Evaluate on Training Set |  60% Train, 40% Test
---|---|---
WordTokenizer |  89.8607 |  57.0222
NGramTokenizer |  89.9043 |  56.5693
Character Unigram |  34.9289 |  33.6113
Character Bigram |  70.5250 |  41.3659
Character Trigram |  94.1979 |  45.6490
TABLE III: The accuracy of SMO classifier with different features

##  V Result

Finally, we used the prepared separate test data set to evaluate our system
via SMO classifier and Character NGramTokenizer filter. Our experimental
result is obtained by using TriGramTokenizer, IDF=False, TF=False. The total
number of the instances is 1350, the accuracy of our model approximately
equals to 39%, and the Mean Absolute Error (MAE) indicating the average value
of the distance between the predicted value of the model and the real value of
the sample is 29%.

Result | Amount | Rate
---|---|---
Total Number of Instances | 1350 | /
Correctly Classified | 531 | 39.35%
Incorrectly Classified | 819 | 60.65%
Mean Absolute Error | 0.2908
TABLE IV: SMO Classified Result

  * •

F1-Score: Precision (PRE) and Recall (REC) are two significant and common
metrics when evaluating the effectiveness of binary classification model.
However, under certain circumstances, Precision and Recall indicators
sometimes affect each other. F1-Measure is generally a harmonic method
balancing and combining both Precision and Recall metrics, which is calculated
as followed:

Class | Precision | Recall | F1-Score | AUC
---|---|---|---|---
MC | 0.426 | 0.625 | 0.506 | 0.716
BE | 0.179 | 0.099 | 0.127 | 0.643
FR | 0.386 | 0.446 | 0.414 | 0.612
CA | 0.440 | 0.345 | 0.386 | 0.706
CH | 0.348 | 0.184 | 0.241 | 0.676
Weighted Avg | 0.377 | 0.393 | 0.373 | 0.673
TABLE V: The values of various indexes of SMO classifier

$\frac{1}{F_{\beta}}=\frac{1}{1+\beta^{2}}(\frac{1}{P}+\frac{\beta^{2}}{P})$
$F_{\beta}=\frac{(\beta^{2}+1)\times P\times R}{(\beta^{2}\times P)+R}$

Where $\beta$ is the parameter, P is Precision and R is Recall. Normally, we
use F1-Score method, which is derived from setting $\beta$ to 1. The range of
F1-Score is from 0 to 1, with 0 represents the worst output. The performance
of the model gets better with the increase of F1-Score.

  * •

ROC Curve: The abscissa of the ROC plane is false positive rate (FPR), and the
ordinate is true positive rate (TPR). For a certain classifier, we can obtain
a TPR and FPR dotted pair according to its performance on the test sample. In
this way, the classifier can be mapped to a point on the ROC plane. By
adjusting the threshold used by this classifier, we can get a curve passing
through (0, 0), (1, 1), which shapes the ROC curve of this classifier.

  * •

AUC:  AUC is a standard metric used to measure the quality of classification
model. The value of AUC is the size of the area below the ROC curve.
Generally, the value of AUC is between 0.5 and 1.0, and a larger AUC value
represents the better performance of the model[4].

Table 5 shows the respective and average accuracy of Precision, Recall,
F1-Score and AUC of five identical classes (i.e.MC, BE, FR, CA, CH) using SMO
classifier.

##  VI Conclusion

Our modified system was built with WEKA data analytic tool and SMO classifier
after measuring and comparing the performance of different classifiers tested
on training set (around 94% accuracy) and 60% for training while the remaining
40% for testing (about 57% accuracy at the best situation). After selecting
SMO as the primary chosen classifier, we then evaluated our system via SMO
classifier and various filter using separate data set. We have figured out
that with our dataset, Character TriGramTokenizer have the best performance.
We found that the reason that we get a lower accuracy is the low distinction
of text and the unselected data, since we used the corpus from Sketch Engine.
Obviously, if the training data set contains mislabeled or appropriated data,
the accuracy of classifier could achieve will be reduced.

## References

  * [1] A. Alshutayri, E. Atwell, A. Alosaimy, J. Dickins, M. Ingleby, and J. Watson, “Arabic language WEKA-based dialect classifier for Arabic automatic speech recognition transcripts,” in Proceedings of the Third Workshop on NLP for Similar Languages, Varieties and Dialects (VarDial3), (Osaka, Japan), pp. 204–211, The COLING 2016 Organizing Committee, Dec. 2016\.
  * [2] S. S. Keerthi, S. K. Shevade, C. Bhattacharyya, and K. R. K. Murthy, “Improvements to platt’s smo algorithm for svm classifier design,” Neural Computation, vol. 13, no. 3, pp. 637–649, 2001.
  * [3] W. Cohen, “Fast effective rule induction,” Twelfth International Conference on Machine Learning: 1995, vol. 95, 10 2000.
  * [4] E. Atwell, C. Brierley, C. Rowland, and J. Anderson, “Semantic pathways: a novel visualisation of varieties of english,” . ICAME Journal of the International Computer Archive of Modern English, vol. 37, pp. 5–36, 01 2013\.