remarkRemark hypothesisHypothesis claimClaim

# Krylov subspace recycling for matrix functions

Liam Burke First author, School Of Mathematics, Trinity College Dublin,
College Green, Dublin 2, Ireland () Email<EMAIL_ADDRESS>Andreas Frommer
Department of Mathematics, Bergische Universität Wuppertal, 42097 Wuppertal,
Germany Gustavo Ramirez-Hidalgo22footnotemark: 2 Kirk M. Soodhalter School Of
Mathematics, Trinity College Dublin, College Green, Dublin 2, Ireland

###### Abstract

We derive an augmented Krylov subspace method with subspace recycling for
computing a sequence of matrix function applications on a set of vectors. The
matrix is either fixed or changes as the sequence progresses. We assume
consecutive matrices are closely related, but make no assumptions on the
relationship between the vectors. We present three versions of the method with
different practical implementations. We demonstrate the effectiveness of the
method using a range of numerical experiments with a selection of functions
and matrices. We primarily focus our attention on the sign function arising in
the overlap formalism of lattice QCD.

###### keywords:

Krylov subspace recycling, matrix functions

65F10, 65F30, 65F50

##  1 Introduction

Many applications in scientific computing require the evaluation of
$f(\textbf{A})\textbf{b}$ where $\textbf{A}\in\mathbb{C}^{n\times n}$ is some
matrix, $\textbf{b}\in\mathbb{C}^{n}$ is a vector, and $f$ is some complex
valued matrix function such that $f(\textbf{A})$ is defined. For example, the
matrix exponential function $f(z)=\mbox{exp}(z)$ arises in the solution of
ordinary differential equations [24, 25]. The matrix logarithm
$f(z)=\mbox{log}(z)$ arises in Markov model analysis [35], and the matrix sign
function $f(z)=\mbox{sign}(z)$ arises in lattice QCD simulations with overlap
fermions [28, 14]. Perhaps the most well known example is the solution of a
linear system where we have $f(z)=\frac{1}{z}$.

In High Performance Computing (HPC) applications the matrix A is large and
sparse, and sometimes not available explicitly. In this setting, Krylov
subspace methods are the only efficient means of computing an approximation to
$f(\textbf{A})\textbf{b}$ by projecting the problem onto a $j$-dimensional
Krylov subspace $\mathcal{V}_{j}$ ($j\ll n$), and treating a smaller
representation of the problem on this subspace via a direct method.

In this paper we propose to treat the action of a matrix function on a vector
via augmented Krylov subspace methods. These methods aid in accelerating the
convergence of standard Krylov subspace methods by augmenting
$\mathcal{V}_{j}$ with an additional $k$ dimensional subspace $\mathcal{U}$
known as the augmentation subspace. The space $\mathcal{U}$ is chosen to
contain useful information to aid in accelerating the convergence of the
method and is commonly taken to be approximations to the $k$ eigenvectors
corresponding to the $k$ eigenvalues closest to a singularity of the function
$f$.

We develop an augmented Krylov subspace approximation to
$f(\textbf{A})\textbf{b}$ by deriving an augmented Krylov subspace
approximation to the shifted linear system appearing in its integral
representation. Although augmented Krylov subspace methods for matrix
functions are not new (see eg. [2, 11]), the purpose of this paper is to
derive a method which allows for subspace recycling between a sequence of $N$
matrix function applications of the form

$f(\textbf{A}^{(i)})\textbf{b}^{(i)},\enspace i=1,2,..,N,$

where $\textbf{A}^{(i)}$ and/or $\textbf{b}^{(i)}$ are changing as the
sequence progresses, and the function $f$ is not restricted to be the inverse
function. Problems of this form arise in simulations of Quantum Chromodynamics
(QCD), the theory of the strong interaction. The overlap formalism of lattice
QCD [31] suffers from computational overhead arising from the above problem
with $f(z)=\mbox{sign}(z)$. Treating this problem by a subspace recycling
method will be more flexible than current methods used in situations where the
matrix remains fixed, and could allow for a more rapid decrease in convergence
of the method as each problem in the sequence is solved. The method in this
paper is the first of its kind proposed to extend subspace recycling to a
sequence of matrix function applications where the matrices in the sequence
are slowly changing.

The layout of this paper is as follows. In Section 2 we give some background
on the literature for augmented Krylov subspace methods for matrix functions.
We also discuss the GCRO-DR algorithm for linear systems [32] which serves as
an inspiration for the method we derive in this work. In Section 3 we discuss
the standard Arnoldi approximation to $f(\textbf{A})\textbf{b}$ and in
particular focus on its derivation using the Cauchy integral representation of
the matrix function. In Section 4 we discuss computational challenges arising
in the overlap formalism of lattice QCD which has served as our main
motivation for this work. In Section 5 we introduce an augmented Krylov
subspace framework for solving the shifted linear systems arising in the
integral representations of $f(\textbf{A})\textbf{b}$. In Section 6 we show
how this framework for shifted systems can be used to derive an augmented
Krylov subspace approximation for $f(\textbf{A})\textbf{b}$. We present three
different versions of a FOM type method, leading to different implementations.
A brief discussion on Stieltjes functions can be found in Section 7. In
Section 8 we derive a harmonic Ritz procedure for our method which differs
slightly from the harmonic Ritz procedure used in GCRO-DR. In Section 9 we
present numerical experiments and test the three versions of the method as
both an augmented method for a single problem, and a recycle method for a
sequence of problems. Brief conclusions are given in Section 10.

##  2 Background

In recent years there has been an interest in improving computational aspects
of the Arnoldi approximation for matrix functions in order to keep up with
application demands and practical computing limitations. A procedure for
implementing restarts into the Arnoldi approximation by computing an error
function using divided differences was first introduced in [10]. In [13] it
was shown that expressing this error function as an integral allows for a more
numerically stable restart procedure based on quadrature. A convergence
analysis was later provided in [12]. A framework for treating the action of a
matrix function on multiple right hand sides using block Krylov subspace
methods was derived in the PhD thesis [29] and can also be found in the
publications [15, 16].

The augmented Krylov subspace methods proposed in the literature aim at
treating problems where the matrix and vector remain fixed. One of the earlier
works has developed a deflation technique which precomputes Schur vectors from
the original matrix and runs a single Arnoldi cycle which orthogonalizes the
Arnoldi vectors against these Schur vectors [2] . Similarly, a procedure for
performing deflation between restarts of the Arnoldi process was introduced in
[11]. This allows one to avoid the computational expense of precomputing a
deflation subspace from the original matrix. The methods in these papers are
based off the Arnoldi-like decomposition [11], and the augmented subspace can
still be treated as a single Krylov subspace. A generalization of the Arnoldi
approximation for non-augmented methods can then be used.

This is similar to ideas presented in the literature for Krylov subspace
methods for the solution of a linear system. Morgan [30] showed that when
GMRES is deflated with harmonic Ritz vectors between restarts, then the
augmented Krylov subspace is still a Krylov subspace, but with a different
starting vector. It thus preserves the usual properties of the Krylov subspace
which allows the approximation to be computed in a manner which looks like the
standard GMRES method. This imposes restrictions on the choice of subspace
$\mathcal{U}$, and in addition, the method can only be used to treat a single
linear system. The GCRO-DR algorithm [32] overcomes these difficulties by
combining ideas from GMRES-DR and GCRO [9] to allow for an arbitrary choice of
$\mathcal{U}$, as well as allowing for recycling between a sequence of
changing linear systems. We refer the reader to [32, 30] for further details
on these methods.

In this paper we aim to extend the ideas in deflated methods for matrix
functions such as [11] to methods which allow for an arbitrary choice of
subspace $\mathcal{U}$ while also allowing for recycling. We will do so using
the main mechanics of the GCRO-DR algorithm for linear systems by maintaining
the relation $\textbf{C}=\textbf{A}\textbf{U}$. We note however, that the
method we develop is not actually based off a GMRES iteration, but rather a
FOM iteration.

##  3 The Arnoldi approximation to $f(\textbf{A})\textbf{b}$

We shortly summarize the essentials on the theoretical justification and the
computation of the well-known Arnoldi approximation to
$f(\textbf{A})\textbf{b}$. We will take inspiration from this approach when
developing the augmented Krylov subspace approximation.

Any matrix function $f(\textbf{A})$ is identical to a matrix polynomial
$p(\textbf{A})$, where $p$ interpolates $f$ at the eigenvalues of A in the
Hermite sense; see [23]. This is why the Krylov subspaces

$\mathcal{K}_{j}(\textbf{A},\textbf{b})=\\{q(\textbf{A})\textbf{b}:\mbox{ $q$
polynomial of degree $<j$}\\}$

are natural spaces from which to choose an appropriate approximation to
$f(\textbf{A})\textbf{b}$. In the case of a linear system,
$f(\textbf{A})=\textbf{A}^{-1}$, a Galerkin condition can be used to define
such an approximation
$\textbf{x}_{j}\in\mathcal{K}_{j}(\textbf{A},\textbf{b})$ via the
orthogonality condition

(1)
$\textbf{b}-\textbf{A}\textbf{x}_{j}\perp\mathcal{K}_{j}(\textbf{A},\textbf{b}).$

In the full orthogonalization method FOM for linear systems [33],
$\textbf{x}_{j}$ from (1) is obtained from the Arnoldi process which computes
a nested orthonormal basis for $\mathcal{K}_{j}(\textbf{A},\textbf{b})$.
Arranging the basis vectors as the columns of a matrix
$\textbf{V}_{j}\in\mathbb{C}^{n\times j}$, the Arnoldi process can be
summarized via the Arnoldi relation

(2)
$\textbf{A}\textbf{V}_{j}=\textbf{V}_{j+1}\overline{\textbf{H}}_{j}=\textbf{V}_{j}\textbf{H}_{j}+h_{j+1,j}\textbf{v}_{j+1}\textbf{e}_{j}^{T},$

where $\overline{\textbf{H}}_{j}\in\mathbb{C}^{(j+1)\times j}$ is an upper
Hessenberg matrix and $\textbf{H}_{j}\in\mathbb{C}^{j\times j}$ arises from
$\overline{\textbf{H}}_{j}$ by deleting the last row. The FOM iterate
$\textbf{x}_{j}$ satisfying (1) can then be computed as
$\textbf{x}_{j}=\|\textbf{b}\|\textbf{V}_{j}\textbf{H}_{j}^{-1}\textbf{e}_{1}$,
where $\textbf{e}_{1}$ is the first canonical unit vector in $\mathbb{C}^{j}$.

A crucial observation for matrix functions is the fact that Krylov subspaces
are invariant under shifts of the matrix by multiples of the identity

$\mathcal{K}_{j}(\textbf{A},\textbf{b})=\mathcal{K}_{j}(\sigma\textbf{I}-\textbf{A},\textbf{b})$

for all $\sigma\in\mathbb{C}$, and that the FOM approximations for
$(\sigma\textbf{I}-\textbf{A})^{-1}\textbf{b}$ are given by

$\textbf{x}_{j}(\sigma)=\|\textbf{b}\|\textbf{V}_{j}(\sigma\textbf{I}-\textbf{H}_{j})^{-1}\textbf{e}_{1},$

with $\textbf{V}_{j}$ and $\textbf{H}_{j}$ built from the Arnoldi process for
A; see e.g. [34].

For functions $f$ which are analytic in an appropriate region, such shifted
systems arise when expressing $f(\textbf{A})$ using Cauchy’s integral formula
as given in the following definition; see [23].

######  Definition 3.1.

For the matrix $\textbf{A}\in\mathbb{C}^{n\times n}$ with spectrum
$\Lambda(\textbf{A})$, if $f$ is analytic in a region containing the closed
contour $\Gamma$ which in turn contains $\Lambda(\textbf{A})$ in its interior,
$f(\textbf{A})$ can be defined as

(3) $f(\textbf{A})=\frac{1}{2\pi
i}\int_{\Gamma}f(\sigma)(\sigma\textbf{I}-\textbf{A})^{-1}d\sigma.$

From this definition we obtain

(4) $f(\textbf{A})\textbf{b}=\frac{1}{2\pi
i}\int_{\Gamma}f(\sigma)(\sigma\textbf{I}-\textbf{A})^{-1}\textbf{b}\hskip
2.84544ptd\sigma,$

which is built from the solutions $\textbf{x}(\sigma)$ of the shifted linear
system

(5) $(\sigma\textbf{I}-\textbf{A})\textbf{x}(\sigma)=\textbf{b}.$

Using the FOM approximations $\textbf{x}_{j}(\sigma)$ for each shift, we then
get the approximation

(6) $f(\textbf{A})\textbf{b}\approx\frac{\|\textbf{b}\|}{2\pi
i}\int_{\Gamma}f(\sigma)\textbf{V}_{j}(\sigma\textbf{I}-\textbf{H}_{j})^{-1}\textbf{e}_{1}\hskip
2.84544ptd\sigma=\|\textbf{b}\|\textbf{V}_{j}f(\textbf{H}_{j})\textbf{e}_{1}.$

Note that the last equality only holds if the spectrum of $\textbf{H}_{j}$
lies inside the contour $\Gamma$.

The last equality is the motivation to take

$\|\textbf{b}\|\textbf{V}_{j}f(\textbf{H}_{j})\textbf{e}_{1}$

as the $j$th Arnoldi approximation to $f(\textbf{A})\textbf{b}$. It is defined
if $f(\textbf{H}_{j})$ is defined, i.e. when $f$ is defined on the spectrum of
$\textbf{H}_{j}$ in the sense of [23]. Note that we can express
$f(\textbf{H}_{j})$ via the Cauchy integral formula using any admissible
contour $\widetilde{\Gamma}$ which encloses the spectrum of $\mathbf{H}_{j}$
but not necessarily that of $\mathbf{A}$.

##  4 The sign function in lattice QCD

This work has been motivated by computational challenges arising in lattice
QCD simulations [27]. The overlap formalism [31] introduces the overlap Dirac
operator

$\textbf{D}_{ovl}=\rho\textbf{I}+\Gamma_{5}\mbox{sign}(\textbf{Q}).$

Here, $\rho$ is a mass parameter, $\textbf{Q}=\Gamma_{5}\textbf{D}$, wherein
the Dirac-Wilson matrix D represents a periodic nearest neighbour coupling on
a 4d space-time lattice with 12 unknowns per lattice site (corresponding to
all possible combinations of four spins and three colors). The matrix
$\Gamma_{5}$ acts identically on the components at each lattice sites by
flipping half of the spins. As a consequence,
$\textbf{Q}=\Gamma_{5}\textbf{D}$ is Hermitian, $\textbf{Q}^{*}=\textbf{Q}$.

The overlap formalism preserves chiral symmetry on the lattice, and this has
significant importance for certain types of observables computed in lattice
QCD simulations. These simulations are thus dependent on overcoming the
algorithmic complexity challenges associated with the overlap discretization
which we will discuss in this section.

Lattice QCD computations with the overlap formalism are dominated by computing
solutions to a sequence of $N$ linear systems of the form

$\textbf{D}_{ovl}^{(i)}\textbf{x}^{(i)}=\textbf{b}^{(i)},\hskip
28.45274pti=1,..,N.$

The computation of quark propagators [27] is a special case wherein the
matrices $\textbf{D}_{ovl}^{(i)}$ remains fixed and the vectors
$\textbf{b}^{(i)}$ are changing from one system to the next. In the initial
phase of generating QCD gauge fields in the Markov chain Monte Carlo process,
the matrix $\textbf{D}_{ovl}^{(i)}$ is slowly changing for every system [6].

The matrices in lattice QCD are large and sparse, and computing an
approximation of the solution for system $i$ can only be done via
(preconditioned) Krylov subspace methods. The Arnoldi algorithm can be used to
build a basis for the Krylov subspace
$\mathcal{K}_{j}(\textbf{D}_{ovl}^{(i)},\textbf{r}_{0}^{(i)})$, where
$\textbf{r}_{0}^{(i)}$ is the initial residual corresponding to an initial
approximation $\textbf{x}_{0}^{(i)}$ for system $i$. Each iteration of the
Arnoldi algorithm requires the application of $\textbf{D}^{(i)}_{ovl}$ to a
new vector. For each system $i$ this requires the evaluation of
$\mbox{sign}(\textbf{Q}^{(i)})\textbf{v}^{(\ell)}$ for a fixed matrix
$\textbf{Q}^{(i)}$ across a sequence of vectors $\textbf{v}^{(\ell)}$,
$\ell=1,..,j$, where $j$ is the Arnoldi cycle length. After each system, the
matrix $\textbf{Q}^{(i)}$ either remains fixed, or changes to a new matrix
$\textbf{Q}^{(i+1)}$, and the process repeats.

The application of the sign function of each matrix on the set of vectors is a
major computational bottleneck in the overlap formalism, and requires an
additional inner Krylov subspace method to compute
$\mbox{sign}(\textbf{Q}^{(i)})$ within each single iterative step in the
solution process for a system with the overlap operator. This makes the
overlap formalism very expensive to work with, and has thus attracted much
interest in developing new methods to address these computational challenges.
We refer the reader to works such as [14], as well as the series of
publications [1, 6, 7, 40]. We also refer the reader to [4] where a Multigrid
preconditioning for the overlap operator was considered. In previous work
[40], exact deflation was used as a means to accelerate the inner iteration.
This requires one to compute a given number of smallest eigenpairs of each
$\textbf{Q}^{(i)}$ to high accuracy before starting the computation of the
sign function. The work in this paper proposes an alternative using Krylov
subspace recycling instead. This does not require one to compute exact small
eigenpairs in a separate process, but interweaves the computation of an
increasingly accurate eigenmodes with the computation of the sign function
itself.

##  5 A derivation of an augmented Krylov subspace approximation to shifted
linear systems

We derive an augmented Krylov subspace approximation for
$f(\textbf{A})\textbf{b}$ based on an integral representation in a similar
manner as for the Arnoldi approximation presented in Section 3 for the non-
augmented case. This involves deriving a common augmented Krylov subspace
approximation method for the family of shifted linear systems (5) appearing in
the integral representation (4). We make use of a residual projection
framework which describes almost all augmented Krylov subspace methods for
standard non-shifted linear systems. The framework was originally proposed in
terms of typical residual constraints in works such as [17, 18, 19, 20], but a
recent survey [37] extends the framework to the case of arbitrary search and
constraint spaces. The work in this section is based off the most general form
of the framework as presented in this survey. The framework describes
augmented Krylov subspace methods broadly in two main steps. The first applies
a standard Krylov subspace method to an appropriately projected problem, and
the second performs an additional projection into the augmentation subspace.
We present an extension of this framework to the case of a family of shifted
linear systems (5). We assume that the initial approximations for all shifts
are zero and that the approximation $\textbf{x}_{j}(\sigma)$ for every shift
$\sigma$ is of the form

$\textbf{x}_{j}(\sigma)=\textbf{s}_{j}(\sigma)+\textbf{t}_{j}(\sigma)$

with $\textbf{t}_{j}(\sigma)$ taken from the same $j$-dimensional Krylov
subspace $\mathcal{V}_{j}$ for all $\sigma$ and $\textbf{s}_{j}(\sigma)$ from
a fixed $k$\- dimensional augmentation subspace $\mathcal{U}$ which does not
depend on the iteration step $j$. We use the $j$ subscript on s since, as we
will show, the contribution from the augmentation subspace will also depend on
$j$. With the columns of U spanning $\mathcal{U}$ and those of
$\textbf{V}_{j}$ spanning $\mathcal{V}_{j}$ the approximation
$\textbf{x}_{j}(\sigma)$ at step $j$ can be expressed as

$\textbf{x}_{j}(\sigma)=\textbf{s}_{j}(\sigma)+\textbf{t}_{j}(\sigma)=\textbf{U}\textbf{z}_{j}(\sigma)+\textbf{V}_{j}\textbf{y}_{j}(\sigma)$

with $\textbf{z}_{j}(\sigma)\in\mathbb{C}^{k}$ and
$\textbf{y}_{j}(\sigma)\in\mathbb{C}^{j}$ to be determined. This is done by
imposing that the residual be orthogonal to the sum of some other
$j$-dimensional Krylov subspace $\widetilde{\mathcal{V}}_{j}$ and some other
augmentation subspace $\widetilde{\mathcal{U}}$ for all shifts $\sigma$,

$\textbf{r}_{j}(\sigma)=\textbf{b}-(\sigma\textbf{I}-\textbf{A})(\textbf{U}\textbf{z}_{j}(\sigma)+\textbf{V}_{j}\textbf{y}_{j}(\sigma))\perp\widetilde{\mathcal{U}}+\widetilde{\mathcal{V}}_{j}.$

This yields the following linear system for $\textbf{z}_{j}(\sigma)$ and
$\textbf{y}_{j}(\sigma)$:

$\begin{bmatrix}\widetilde{\textbf{U}}^{*}(\sigma\textbf{I}-\textbf{A})\textbf{U}&\widetilde{\textbf{U}}^{*}(\sigma\textbf{I}-\textbf{A})\textbf{V}_{j}\\\
\widetilde{\textbf{V}}^{*}_{j}(\sigma\textbf{I}-\textbf{A})\textbf{U}&\widetilde{\textbf{V}}^{*}_{j}(\sigma\textbf{I}-\textbf{A})\textbf{V}_{j}\end{bmatrix}\begin{bmatrix}\textbf{z}_{j}(\sigma)\\\
\textbf{y}_{j}(\sigma)\end{bmatrix}=\begin{bmatrix}\widetilde{\textbf{U}}^{*}\textbf{b}\\\
\widetilde{\textbf{V}}_{j}^{*}\textbf{b}\end{bmatrix}.$

Eliminating the bottom left block we can express $\textbf{y}_{j}(\sigma)$ as
the solution of

(7)
$\widetilde{\textbf{V}}_{j}^{*}(\textbf{I}-\textbf{Q}_{\sigma})(\sigma\textbf{I}-\textbf{A})\textbf{V}_{j}\textbf{y}_{j}(\sigma)=\widetilde{\textbf{V}}_{j}^{*}(\textbf{I}-\textbf{Q}_{\sigma})\textbf{b},$

with the projector
$\textbf{Q}_{\sigma}:=(\sigma\textbf{I}-\textbf{A})\textbf{U}(\widetilde{\textbf{U}}^{*}(\sigma\textbf{I}-\textbf{A})\textbf{U})^{-1}\widetilde{\textbf{U}}^{*}$.
The coefficients $\textbf{z}_{j}(\sigma)$ for the contribution from the
augmenting subspace are then obtained from $\textbf{y}_{j}(\sigma)$ as

(8)
$\textbf{z}_{j}(\sigma)=(\widetilde{\textbf{U}}^{*}(\sigma\textbf{I}-\textbf{A})\textbf{U})^{-1}\widetilde{\textbf{U}}^{*}\textbf{b}-(\widetilde{\textbf{U}}^{*}(\sigma\textbf{I}-\textbf{A})\textbf{U})^{-1}\widetilde{\textbf{U}}^{*}(\sigma\textbf{I}-\textbf{A})\textbf{V}_{j}\textbf{y}_{j}(\sigma),$

From (7) we see that computing $\textbf{y}_{j}(\sigma)$ is equivalent to
applying the following projection method.

Find $\textbf{t}_{j}(\sigma)\in\mathcal{V}_{j}$ as an approximate solution
corresponding to shift $\sigma$ for the projected and shifted linear system
(9)
$(\textbf{I}-\textbf{Q}_{\sigma})(\sigma\textbf{I}-\textbf{A})\textbf{x}(\sigma)=(\textbf{I}-\textbf{Q}_{\sigma})\textbf{b}$
such that
$\textbf{r}_{j}(\sigma)=(\textbf{I}-\textbf{Q}_{\sigma})(\textbf{b}-(\sigma\textbf{I}-\textbf{A})\textbf{t}_{j}(\sigma))\perp\widetilde{\mathcal{V}}_{j}$.

As outlined in the framework [18, 37] a natural first choice one may consider
for $\mathcal{V}_{j}$ is then the Krylov subspace corresponding to the
projected problem

$\mathcal{K}_{j}((\textbf{I}-\textbf{Q}_{\sigma})(\sigma\textbf{I}-\textbf{A}),(\textbf{I}-\textbf{Q}_{\sigma})\textbf{b}).$

However this Krylov subspace is built from a projected operator which is a
function of the shift $\sigma$, and thus shift invariance no longer holds, and
a practical implementation would require one to build a separate basis for
every shift. This is one of the main difficulties one encounters when
attempting to develop a subspace recycling method for the solution of shifted
systems [36, 38].

One way to overcome this difficulty is to use the observation made in [5]
where it was shown that the projected problem (9) has an unprojected
equivalent given by

Find $\textbf{t}_{j}(\sigma)\in\mathcal{V}_{j}$ as an approximate solution
corresponding to shift $\sigma$ for the shifted linear system
$(\sigma\textbf{I}-\textbf{A})\textbf{x}(\sigma)=\textbf{b}$ such that
$\textbf{r}_{j}(\sigma)=\textbf{b}-(\sigma\textbf{I}-\textbf{A})\textbf{t}_{j}(\sigma)\perp(\textbf{I}-\textbf{Q}_{\sigma})^{*}\widetilde{\mathcal{V}}_{j}$.

The proof of this equivalence is a trivial extension of the proof shown in [5]
to the shifted case and is thus not repeated here. We now propose to solve the
unprojected problem using the unprojected, shift-invariant Krylov subspace
$\mathcal{K}_{j}(\textbf{A},\textbf{b})$.

The projector $(\textbf{I}-\textbf{Q}_{\sigma})$ is now built into the
constraint space and only appears in the equation one must solve for the
coefficients $\textbf{y}_{j}(\sigma)$ given by (7). Once we have computed
$\textbf{y}_{j}(\sigma)$ we can incorporate $\textbf{z}_{j}(\sigma)$ from (8)
implicitly into the full solution approximation to get

(10)
$\textbf{x}_{j}(\sigma)=\textbf{V}_{j}\textbf{y}_{j}(\sigma)+\textbf{U}(\widetilde{\textbf{U}}^{*}(\sigma\textbf{I}-\textbf{A})\textbf{U})^{-1}\widetilde{\textbf{U}}^{*}(\textbf{b}-(\sigma\textbf{I}-\textbf{A})\textbf{V}_{j}\textbf{y}_{j}(\sigma)).$

This approximation serves a general augmented Krylov subspace approximation
for shifted systems of the form (5) which fits into the framework in [18, 37]
by splitting the problem for computing $\textbf{y}_{j}(\sigma)$ and
$\textbf{z}_{j}(\sigma)$ into two separate problems. Alternatively, one can
treat $\textbf{y}_{j}(\sigma)$ and $\textbf{z}_{j}(\sigma)$ as a single
unknown $\overline{\textbf{y}}_{j}(\sigma)$ defined by

$\displaystyle\overline{\textbf{y}}_{j}(\sigma):=\begin{bmatrix}\textbf{z}_{j}(\sigma)\\\
\textbf{y}_{j}(\sigma)\end{bmatrix}\in\mathbb{C}^{k+j}.$

It can then be shown that the solution approximation (10) is equivalent to

(11)
$\textbf{x}_{j}(\sigma)=\widehat{\textbf{V}}_{j}\begin{bmatrix}\widetilde{\textbf{W}}^{*}_{j}(\sigma\textbf{I}-\textbf{A})\widehat{\textbf{V}}_{j}\end{bmatrix}^{-1}\widetilde{\textbf{W}}^{*}_{j}\textbf{b}$

with
$\widehat{\textbf{V}}_{j}=\begin{bmatrix}\textbf{U}&\textbf{V}_{j}\end{bmatrix}\in\mathbb{C}^{n\times(k+j)}$
and
$\widetilde{\textbf{W}}_{j}=\begin{bmatrix}\widetilde{\textbf{U}}&\widetilde{\textbf{V}}_{j}\end{bmatrix}\in\mathbb{C}^{n\times(k+j)}$.

This gives two mathematically equivalent means of computing the same
approximation $\textbf{x}_{j}(\sigma)$ to the solution of the shifted system
(5) in terms of arbitrary search and constraint spaces. The first one (10) was
derived directly from the framework in [18, 37], while the second (11) treated
the two corrections as a single unknown to be computed at once.

One can use this general framework to derive any new augmented unprojected
Krylov subspace approximation for solving families of shifted system (5) using
appropriate choices of the subspaces $\widetilde{\mathcal{V}}_{j}$,
$\mathcal{U}$ and $\widetilde{\mathcal{U}}$ while keeping
$\mathcal{V}_{j}=\mathcal{K}_{j}(\textbf{A},\textbf{b})$. The choice
$\widetilde{\mathcal{V}}_{j}=\textbf{A}\mathcal{K}(\textbf{A},\textbf{b})$ and
$\widetilde{\mathcal{U}}=\textbf{A}\mathcal{U}$ leads to an unprojected
augmented GMRES-type method for shifted systems, while the choice
$\widetilde{\mathcal{V}}_{j}=\mathcal{K}(\textbf{A},\textbf{b})$ and
$\widetilde{\mathcal{U}}=\mathcal{U}$ leads to an unprojected augmented FOM-
type method which can be viewed as an extension of the unprojected augmented
FOM method derived in [5] to the case of shifted systems of the form (5).
Since we are dealing with matrix functions, where the standard Arnoldi
approximation is of FOM-type, we now focus on FOM type methods and make these
choices.

##  6 Recycled FOM for functions of matrices (rFOM2)

Taking the results from the previous section and using the integral
representation (4) we explicitly define the augmented unprojected FOM
approximation of $f(\textbf{A})\textbf{b}$ with subspaces
$\widetilde{\mathcal{V}}_{j}=\mathcal{V}_{j}=\mathcal{K}_{j}(\textbf{A},\textbf{b})$
and $\widetilde{\mathcal{U}}=\mathcal{U}$ as

(12) $\displaystyle f(\textbf{A})\textbf{b}\approx\frac{1}{2\pi
i}\int_{\Gamma}f(\sigma)\textbf{x}_{j}(\sigma)\hskip 2.84544ptd\sigma,$

where $\textbf{x}_{j}(\sigma)$ is the shifted augmented FOM approximation to
the solution of the shifted system
$(\sigma\textbf{I}-\textbf{A})\textbf{x}=\textbf{b}$ given in (10) or,
equivalently, (11). The use case that we have in mind is that we must compute
a sequence of expressions of the form $f(\textbf{A})\textbf{b}$ with varying b
and/or A and that we use information acquired from the computation of one
expression to obtain good choices for the augmenting space $\mathcal{U}$. This
is why, keeping with the notation of [15, 16, 29] we denote our method as
rFOM2 (recycled FOM for Functions Of Matrices).

Depending on how we treat (12) computationally, and in particular on how we
use quadrature rules on the different integrals occurring in equivalent
representations of (12), we obtain three versions of the FOM-type method that
we discuss below. We use the notation $n_{quad}$ to denote the number of
quadrature points used, and $z_{\ell}$ and $\omega_{\ell}$,
$\ell=1,\ldots,n_{\textrm{quad}}$ to denote the nodes and the weights,
respectively, of a quadrature rule for the contour $\Gamma$ which yields the
approximation

$\sum_{\ell=1}^{n_{\textrm{quad}}}\omega_{\ell}f(z_{\ell})\textbf{x}_{j}(z_{\ell})\approx\frac{1}{2\pi
i}\int_{\Gamma}f(\sigma)\textbf{x}_{j}(\sigma)\hskip 2.84544ptd\sigma.$

Approximating $f(\textbf{A})\textbf{b}$ in its integral form via quadrature
has been discussed in works such as [3, 21]. In addition, numerical quadrature
also plays a role in numerically stable restart procedures in the Arnoldi
approximation [13]. We stress again that our three versions differ because
they use numerical integration for different integrals.

Due to the nature of the expressions we derive inside the Cauchy integral, it
is not possible to modify the contour to just contain the spectrum of a
smaller matrix as with the standard Arnoldi approximation; see Section 3. As
we see it may be necessary to enlarge the contour to enclose additional points
depending on the choice of $\mathcal{U}$. This is a potential limitation of
our method worth considering, since it requires a priori information on the
spectrum of each matrix in order to do the quadrature approximation. We note
that in our numerical experiments we constructed a contour which we knew only
enclosed the spectrum of $\textbf{H}_{j}$, and this did not pose any issues.

The common feature of all three versions is that we perform the Arnoldi
process to compute the orthogonal basis of
$\mathcal{V}_{j}=\mathcal{K}(\textbf{A},\textbf{b})$. This requires $j$
matrix-vector multiplications and $\mathcal{O}(nj^{2})$ arithmetic work for
the orthogonalization. The three versions differ in the integrals that are
evaluated numerically, and we discuss the cost of each version below. For
further implementation details, we refer the reader to the MATLAB code111Code
available at https://github.com/burkel8/MatrixFunctionRecycling.

###  6.1 rFOM2 Version 1

The first version of rFOM2 comes from the framework presented in [18, 37]. We
apply the quadrature rule directly to (12) using the expression (10) for
$\textbf{x}_{j}(\sigma)$ (with $\widetilde{\mathcal{U}}=\mathcal{U}$ and
$\widetilde{\mathcal{V}}_{j}=\mathcal{K}_{j}(\textbf{A},\textbf{b})$).
Expression (10) first requires solving equation (7) for
$\textbf{y}_{j}(\sigma)$ where the projector $\textbf{Q}_{\sigma}$ is now

$\textbf{Q}_{\sigma}=(\sigma\textbf{U}-\textbf{C})(\sigma\textbf{U}^{*}\textbf{U}-\textbf{U}^{*}\textbf{C})^{-1}\textbf{U}^{*}.$

This can be done efficiently using the Arnoldi relation which yields

(13)
$(\textbf{V}_{j}^{*}(\textbf{I}-\textbf{Q}_{\sigma})\textbf{V}_{j+1}(\sigma\overline{\textbf{I}}-\overline{\textbf{H}}_{j}))\textbf{y}_{j}(\sigma)=\textbf{V}^{*}_{j}(\textbf{I}-\textbf{Q}_{\sigma})\textbf{b}.$

Here $\overline{\textbf{I}}\in\mathbb{R}^{(j+1)\times j}$ is a $j\times j$
identity matrix padded with an additional row of zeros. The expression (10) is
then evaluated as

$\textbf{x}_{j}(\sigma)=\textbf{V}_{j}\textbf{y}_{j}(\sigma)+\textbf{U}(\sigma\textbf{U}^{*}\textbf{U}-\textbf{U}^{*}\textbf{C})^{-1}(\textbf{U}^{*}\textbf{b}-\textbf{U}^{*}\textbf{V}_{j+1}(\sigma\overline{\textbf{I}}-\overline{\textbf{H}}_{j})\textbf{y}_{j}(\sigma)).$

With these choices, our first version of rFOM2 thus computes the approximation
$\widetilde{f}_{1}$ to $f(\textbf{A})\textbf{b}$ as

$\displaystyle\widetilde{f}_{1}=$
$\displaystyle\textbf{V}_{j}\sum_{\ell=1}^{n_{\textrm{quad}}}\omega_{\ell}f(z_{\ell})\textbf{y}_{j}(z_{\ell})$
$\displaystyle+\textbf{U}\sum_{\ell=1}^{n_{\textrm{quad}}}\omega_{\ell}f(z_{\ell})(z_{\ell}\textbf{U}^{*}\textbf{U}-\textbf{U}^{*}\textbf{C})^{-1}(\textbf{U}^{*}\textbf{b}-\textbf{U}^{*}\textbf{V}_{j+1}(z_{\ell}\overline{\textbf{I}}-\overline{\textbf{H}}_{j})\textbf{y}_{j}(z_{\ell})).$

For computational efficiency, it is appropriate to compute the coefficient
matrix for $\textbf{y}_{j}$ in (13) as

(14)
$\textbf{V}_{j}^{*}(\textbf{I}-\textbf{Q}_{\sigma})\textbf{V}_{j+1}(\sigma\overline{\textbf{I}}-\overline{\textbf{H}}_{j})=\sigma\textbf{I}-\textbf{H}_{j}-\textbf{K}(\sigma)\textbf{L}(\sigma)\textbf{M}(\sigma\overline{\textbf{I}}-\overline{\textbf{H}}_{j})\
$

with

$\textbf{K}(\sigma)=\sigma\textbf{V}_{j}^{*}\textbf{U}-\textbf{V}_{j}^{*}\textbf{C},\enspace\textbf{L}(\sigma)=(\sigma\textbf{U}^{*}\textbf{U}-\textbf{U}^{*}\textbf{C})^{-1},\enspace\textbf{M}=\textbf{U}^{*}\textbf{V}_{j+1},$

which allows to avoid redundant computations.

This algorithm requires non-singularity of the matrix
$\textbf{U}^{*}(\sigma\textbf{I}-\textbf{A})\textbf{U}$. This can be
guaranteed provided that our contour of integration encloses the field of
values
$\mathcal{F}(\textbf{A}):=\\{\textbf{t}^{*}\textbf{A}\textbf{t},\textbf{t}\in\mathbb{C},\|\textbf{t}\|_{2}=1\\}$
of A which is a convex superset of the spectrum of A; see [26] . We have the
following proposition.

######  Proposition 6.1.

If U is of full rank and the contour $\Gamma$ encloses
$\mathcal{F}(\textbf{A})$, then the matrix
$\textbf{U}^{*}(\sigma\textbf{I}-\textbf{A})\textbf{U}$ is non-singular for
all $\sigma\in\Gamma$.

######  Proof 6.2.

Assume there is some non-zero $\textbf{p}\in\mathbb{C}^{k}$ such that
$\textbf{U}^{*}(\sigma\textbf{I}-\textbf{A})\textbf{U}\textbf{p}=\textbf{0}$.
Let $\textbf{U}=\textbf{Z}\textbf{S}$ be the QR-factorization of U where
$\textbf{Z}\in\mathbb{C}^{n\times k}$ has orthonormal columns and
$\textbf{S}\in\mathbb{C}^{k\times k}$ is non-singular since U has full rank.
The equality
$\textbf{U}^{*}(\sigma\textbf{I}-\textbf{A})\textbf{U}\textbf{p}=\textbf{0}$
implies
$\textbf{p}^{*}\textbf{U}^{*}(\sigma\textbf{I}-\textbf{A})\textbf{U}\textbf{p}=\textbf{0}$
and thus
$\sigma-\textbf{q}^{*}\textbf{Z}^{*}\textbf{A}\textbf{Z}\textbf{q}=\textbf{0}$
for the unit vector Zq with
$\textbf{q}=\tfrac{1}{\|\textbf{S}\textbf{p}\|}\textbf{S}\textbf{p}$. But this
would mean that $\sigma\in\mathcal{F}(\textbf{A})$, which is impossible for
$\sigma\in\Gamma.$

Algorithm 1 rFOM2 \- A recycled FOM approximation for
$f(\textbf{A})\textbf{b}$ (Version 1)

1: Input: $\textbf{A}\in\mathbb{C}^{n\times n}$,
$\textbf{b}\in\mathbb{C}^{n}$, scalar function $f$,
$\textbf{U}\in\mathbb{C}^{n\times k}$ whose columns span the augmenting
recycled subspace $\mathcal{U}$, $\textbf{C}=\textbf{A}\textbf{U}$, Arnoldi
cycle length $j$

2: Build a basis for $\mathcal{K}_{j}(\textbf{A},\textbf{b})$ via the Arnoldi
process generating $\textbf{V}_{j+1}$ and $\overline{\textbf{H}}_{j}$

3: Determine contour containing the spectrum of A and determine quadrature
rule with nodes $z_{\ell}$ and weights
$\omega_{\ell},\ell=1,\ldots,n_{\textrm{quad}}$

4: Set $\textbf{t}_{1}=0\in\mathbb{C}^{j}$,
$\textbf{t}_{2}=0\in\mathbb{C}^{k}$

5: for $\ell=1,..,n_{\textrm{quad}}$ do

6: $\mu\leftarrow\omega_{\ell}f(z_{\ell})$

7: Solve
$(z_{\ell}\textbf{I}-\textbf{H}_{j}-\textbf{K}(z_{\ell})\textbf{L}(z_{\ell})\textbf{M}(z_{\ell}\overline{\textbf{I}}-\overline{\textbf{H}}_{j}))\textbf{y}=\textbf{V}^{*}_{j}(\textbf{I}-\textbf{Q})\textbf{b}$
for y

8: $\textbf{t}_{1}\leftarrow\textbf{t}_{1}+\mu\hskip 2.84544pt\textbf{y}$

9: $\textbf{t}_{2}\leftarrow\textbf{t}_{2}+\mu\hskip
2.84544pt(z_{\ell}\textbf{U}^{*}\textbf{U}-\textbf{U}^{*}\textbf{C})^{-1}(\textbf{U}^{*}\textbf{b}-\textbf{U}^{*}\textbf{V}_{j+1}(z_{\ell}\overline{\textbf{I}}-\overline{\textbf{H}}_{j})\textbf{y})$

10: end for

11:
$\widetilde{f}_{1}\leftarrow\textbf{V}_{j}\textbf{t}_{1}+\textbf{U}\textbf{t}_{2}$

12: Update U and C such that $\textbf{C}=\textbf{A}\textbf{U}$ to use for next
matrix A

Table 1 gives details on the arithmetic cost of an efficient implementation of
the algorithm, excluding the matrix-vector products with A. It reports the
highest order leading term in the respective arithmetic work in units of one
multiplication plus one addition. We do not count matrix scalings, matrix
additions and matrix-vector products since their cost is of lower order than
that for matrix products and inversions. We treat the Hessenberg matrices as
full matrices, and we assume that we evaluate the matrix products in (14) from
right to left, and that we compute the inverse $\textbf{L}(\sigma)$ explicitly
(with work $k^{3}$). Note that we don’t count any work for line 9, since the
computations there can be arranged as matrix-vector products (and vector
additions) only.

operation | arithmetic cost | how often
---|---|---
$\textbf{U}^{*}\textbf{U}$ and $\textbf{U}^{*}\textbf{C}$ | $2k^{2}n$ | 1
$\textbf{V}_{j}^{*}\textbf{U}$ and $\textbf{V}_{j}^{*}\textbf{C}$ | $2kjn$ | 1
$\textbf{M}=\textbf{U}^{*}\textbf{V}_{j+1}$ | $k(j+1)n$ | 1
$\textbf{M}(z_{\ell}\textbf{I}-\overline{\textbf{H}}_{j}):=\textbf{R}(z_{\ell})$ | $kj(j+1)$ | $n_{quad}$
$\textbf{L}(z_{\ell})=(z_{\ell}\textbf{U}^{*}\textbf{U}-\textbf{U}^{*}\textbf{C})^{-1}$ | $k^{3}$ | $n_{\textrm{quad}}$
$\textbf{L}(z_{\ell})\textbf{R}(z_{\ell})=:\textbf{S}(z_{\ell})$ | $k^{2}j$ | $n_{quad}$
$\textbf{K}(z_{\ell})\textbf{S}(z_{\ell})=:\textbf{T}(z_{\ell})$ | $kj^{2}$ | $n_{quad}$
Solve for y, system matrix is $\textbf{T}(z_{\ell})$ | $\frac{1}{3}j^{3}$ | $n_{quad}$
total cost:
$n(2k^{2}+2jk+k(j+1))+n_{quad}(kj(j+1)+k^{3}+k^{2}j+kj^{2}+\frac{1}{3}j^{3})$
Table 1: Arithmetic cost of rFOM2 version 1.

###  6.2 rFOM2 Version 2

We present Algorithm 2 as a potential alternative to Algorithm 1. It is
derived from equation (11) and treats the contribution from the Krylov
subspace and the augmentation space as a single unknown. Version 2 makes use
of an augmented Arnoldi relation as shown in the following proposition. This
relation contains a scaling matrix $\textbf{D}\in\mathbb{C}^{k\times k}$ which
has been included for generality. Its main use, however, is in Algorithm 3 and
will be discussed later in the text.

######  Proposition 6.3.

For the choices of subspaces
$\widetilde{\mathcal{V}}_{j}=\mathcal{V}_{j}=\mathcal{K}_{j}(\textbf{A},\textbf{b})$
and $\widetilde{\mathcal{U}}=\mathcal{U}$, assume we then scale U via
$\textbf{U}\leftarrow\textbf{U}\textbf{D}$ where D is a $k\times k$ diagonal
matrix. Then the augmented Krylov subspace approximation (11) to the shifted
system $(\sigma\textbf{I}-\textbf{A})\textbf{x}=\textbf{b}$ can be written as

(15)
$\textbf{x}_{j}(\sigma)=\widehat{\textbf{V}}_{j}(\widehat{\textbf{V}}_{j}^{*}\widehat{\textbf{W}}_{j}(\sigma\textbf{I}-\textbf{G}_{j})+\widehat{\textbf{V}}^{*}_{j}\textbf{R}_{\sigma})^{-1}\widehat{\textbf{V}}_{j}^{*}\textbf{b},$

with the augmented Hessenberg matrix

(16) $\textbf{G}_{j}=\begin{bmatrix}\textbf{D}&\textbf{0}\\\
\textbf{0}&\textbf{H}_{j}\end{bmatrix}\in\mathbb{C}^{(k+j)\times(k+j)},$

and

$\displaystyle\textbf{R}_{\sigma}=\begin{bmatrix}\sigma(\textbf{U}-\textbf{C})&-h_{j+1,j}\textbf{v}_{j+1}\textbf{e}_{j}^{T}\end{bmatrix}\in\mathbb{C}^{n\times(k+j)},\enspace\widehat{\textbf{W}}_{j}=\begin{bmatrix}\textbf{C}&\textbf{V}_{j}\end{bmatrix},\enspace\textbf{C}=\textbf{A}\textbf{U}.$

######  Proof 6.4.

After the scaling $\textbf{U}\leftarrow\textbf{U}\textbf{D}$ we have

AU $\displaystyle=$ $\displaystyle\textbf{C}\textbf{D},$
$\displaystyle(\sigma\textbf{I}-\textbf{A})\textbf{U}$ $\displaystyle=$
$\displaystyle\textbf{C}(\sigma\textbf{I}-\textbf{D})+\sigma(\textbf{U}-\textbf{C}).$

Combining this with the Arnoldi relation (2) yields the shifted augmented
Arnoldi relation

(17)
$\displaystyle(\sigma\textbf{I}-\textbf{A})\begin{bmatrix}\textbf{U}&\textbf{V}_{j}\end{bmatrix}$
$\displaystyle=$
$\displaystyle\begin{bmatrix}\textbf{C}&\textbf{V}_{j}\end{bmatrix}\begin{bmatrix}\sigma\textbf{I}-\textbf{D}&0\\\
0&\sigma\textbf{I}-\textbf{H}_{j}\end{bmatrix}+\begin{bmatrix}\sigma(\textbf{U}-\textbf{C})&-h_{j+1,j}\textbf{v}_{j+1}\textbf{e}_{j}^{T}\end{bmatrix}.$

The result now follows from (11) by replacing the term
$(\sigma\textbf{I}-\textbf{A})\begin{bmatrix}\textbf{U}&\textbf{V}_{j}\end{bmatrix}$
with the right hand side of (17), and observing that in (11) we also have
$\widehat{\textbf{V}}_{j}=\widetilde{\textbf{W}}_{j}$ since
$\mathcal{U}=\widetilde{\mathcal{U}}$.

Using (15) in the integral representation gives the approximation

(18) $f(\textbf{A})\textbf{b}\approx\widehat{\textbf{V}}_{j}\frac{1}{2\pi
i}\int_{\Gamma}f(z)(\widehat{\textbf{V}}_{j}^{*}\widehat{\textbf{W}}_{j}(\sigma\textbf{I}-\textbf{G}_{j})+\widehat{\textbf{V}}^{*}_{j}\textbf{R}_{\sigma})^{-1}d\sigma\widehat{\textbf{V}}_{j}^{*}\textbf{b},$

which, using a quadrature rule, gives our second version of rFOM2, giving the
approximation

$\widetilde{f}_{2}=\widehat{\textbf{V}}_{j}\sum_{\ell=1}^{n_{\textrm{quad}}}\omega_{\ell}f(z_{\ell})(\widehat{\textbf{V}}_{j}^{*}\widehat{\textbf{W}}_{j}(z_{\ell}\textbf{I}-\textbf{G}_{j})+\widehat{\textbf{V}}^{*}_{j}\textbf{R}_{z_{\ell}})^{-1}\widehat{\textbf{V}}_{j}^{*}\textbf{b}.$

The resulting algorithm is given as Algorithm 2.

Algorithm 2 rFOM2 \- A recycled FOM approximation for
$f(\textbf{A})\textbf{b}$ (Version 2)

1: Input: $\textbf{A}\in\mathbb{C}^{n\times n}$,
$\textbf{b}\in\mathbb{C}^{n}$, scalar function $f$,
$\textbf{U}\in\mathbb{C}^{n\times k}$ whose columns span the augmenting
recycled subspace $\mathcal{U}$, $\textbf{C}=\textbf{A}\textbf{U}$, Arnoldi
cycle length $j$

2: Build a basis for $\mathcal{K}_{j}(\textbf{A},\textbf{b})$ via the Arnoldi
process, generating $\textbf{V}_{j+1}$ and $\overline{\textbf{H}}_{j}$

3: Choose scaling matrix $\textbf{D}\in\mathbb{C}^{k\times k}$

4: $\textbf{U}\leftarrow\textbf{U}\textbf{D}$,
$\widehat{\textbf{V}}_{j}\leftarrow\begin{bmatrix}\textbf{U}&\textbf{V}_{j}\end{bmatrix}$,
$\widehat{\textbf{W}}_{j}\leftarrow\begin{bmatrix}\textbf{C}&\textbf{V}_{j}\end{bmatrix}$

5: Determine contour containing the spectrum of A and determine quadrature
rule for the integral in (18) with nodes $z_{\ell}$ and weights
$\omega_{\ell},\ell=1,\ldots,n_{\textrm{quad}}$

6: Set $\textbf{t}=0\in\mathbb{C}^{j+k}$

7: for $\ell=1,..,n_{\textrm{quad}}$ do

8: $\mu\leftarrow\omega_{\ell}f(z_{\ell})$

9:
$\textbf{t}\leftarrow\textbf{t}+\mu\big{(}\widehat{\textbf{V}}_{j}^{*}\widehat{\textbf{W}}_{j}(z_{\ell}\textbf{I}-\textbf{G}_{j})+\widehat{\textbf{V}}^{*}_{j}\textbf{R}_{z_{\ell}}\big{)}^{-1}\widehat{\textbf{V}}^{*}_{j}\textbf{b}$

10: end for

11: $\widetilde{f}_{2}\leftarrow\widehat{\textbf{V}}_{j}\textbf{t}$

12: Update U and C such that $\textbf{C}=\textbf{A}\textbf{U}$ to use for next
matrix A and vector b

The attractive feature of version 2 is that for each quadrature point only one
inversion of a $(j+k)\times(j+k)$ matrix is needed. A cost summary of version
$2$ is given in Table 2. We use the same approach as for the summary in Table
1 and count only matrix-matrix operations. We count the cost of constructing
$\widehat{\textbf{V}}_{j}^{*}\textbf{R}_{z_{\ell}}$ in terms of the cost of
constructing its block components

$\widehat{\textbf{V}}^{*}_{j}\textbf{R}_{z_{\ell}}=\begin{bmatrix}z_{\ell}\textbf{U}^{*}(\textbf{U}-\textbf{C})&\textbf{U}^{*}\textbf{M}\\\
z_{\ell}\textbf{V}_{j}^{*}(\textbf{U}-\textbf{C})&\textbf{V}_{j}^{*}\textbf{M}\end{bmatrix}$

where $\textbf{M}=-h_{j+1,j}\textbf{v}_{j+1}\textbf{e}_{j}^{T}$ is a rank 1
matrix. This means the products $\textbf{U}^{*}\textbf{M}$ and
$\textbf{V}^{*}_{j}\textbf{M}$ are obtained as matrix-vector and vector-vector
operations and thus do not contribute to leading order terms in the arithmetic
cost.

operation | arithmetic cost | how often
---|---|---
$\widehat{\textbf{V}}_{j}^{*}\widehat{\textbf{W}}_{j}$ | $n(j+k)^{2}$ | 1
$\textbf{U}^{*}(\textbf{U}-\textbf{C})$ | $nk^{2}$ | 1
$\textbf{V}_{j}^{*}(\textbf{U}-\textbf{C})$ | $nkj$ | 1
$(\widehat{\textbf{V}}_{j}^{*}\widehat{\textbf{W}}_{j})(z_{\ell}\textbf{I}-\textbf{G}_{j})$ | $(j+k)^{3}$ | $n_{quad}$
Linear solve with $(\widehat{\textbf{V}}_{j}^{*}\widehat{\textbf{W}}_{j}(z_{\ell}\textbf{I}-\textbf{G}_{j})+\widehat{\textbf{V}}_{j}^{*}\textbf{R}_{z_{\ell}})$ | $\frac{1}{3}(j+k)^{3}$ | $n_{quad}$
total cost: $n((j+k)^{2}+k^{2}+kj)+\frac{4}{3}n_{quad}(j+k)^{3}$
Table 2: Arithmetic cost of rFOM2 version 2.

###  6.3 rFOM2 Version 3

Our third version of rFOM2 is based on a reformulation of (18) which involves
an evaluation of $f$ on the extended Hessenberg matrix $\textbf{G}_{j}$ from
(16). Approximation in terms of an evaluation of $f$ on some small Hessenberg
matrix is typical of how standard Krylov subspace methods for matrix functions
are implemented. We have seen this as the case for the plain Arnoldi
approximation, and it is also the case for the deflated Arnoldi like
approximation in [11]. As we will see in the following proposition, in our
situation we cannot write the full approximation completely in terms of
$f(\textbf{G}_{j})$. We can however, write it as the sum of a term involving
$f(\textbf{G}_{j})$ and another integral expression. The purpose of this
formulation is to introduce a term which can be evaluated without numerical
quadrature. This takes away much of the dependence of the approximation on the
quadrature and allows us to achieve a desired level of accuracy with less
quadrature points as shown in the numerical experiments in Section 9. We
derive version 3 of rFOM2 using the following proposition.

######  Proposition 6.5.

Assume that the contour $\Gamma$ not only encloses the spectrum of A, but also
that of $\textbf{G}_{j}$, then the approximation (18) can equivalently be
written as

(19)
$f(\textbf{A})\textbf{b}\approx\widehat{\textbf{V}}_{j}f(\textbf{G}_{j})(\widehat{\textbf{V}}_{j}^{*}\widehat{\textbf{W}}_{j})^{-1}\widehat{\textbf{V}}_{j}^{*}\textbf{b}-\widehat{\textbf{V}}_{j}\hskip
0.28436pt\mathcal{I}\hskip 0.28436pt\widehat{\textbf{V}}_{j}^{*}\textbf{b}$

with

$\displaystyle\mathcal{I}$ $\displaystyle=$ $\displaystyle\frac{1}{2\pi
i}\int_{\Gamma}f(\sigma)(\widehat{\textbf{V}}^{*}_{j}\widehat{\textbf{W}}_{j}(\sigma\textbf{I}-\textbf{G}_{j}))^{-1}\textbf{S}_{j}(\sigma)\widehat{\textbf{V}}_{j}^{*}\textbf{R}_{\sigma}(\widehat{\textbf{V}}^{*}_{j}\widehat{\textbf{W}}_{j}(\sigma\textbf{I}-\textbf{G}_{j}))^{-1}\hskip
2.84544ptd\sigma,$

where
$\textbf{S}_{j}(\sigma)\;=\;(\textbf{I}+\widehat{\textbf{V}}_{j}^{*}\textbf{R}_{\sigma}(\widehat{\textbf{V}}^{*}_{j}\widehat{\textbf{W}}_{j}(\sigma\textbf{I}-\textbf{G}_{j}))^{-1})^{-1}$.

######  Proof 6.6.

By the Sherman–Morrison–Woodbury identity we have

$\displaystyle(\widehat{\textbf{V}}_{j}^{*}\widehat{\textbf{W}}_{j}(\sigma\textbf{I}-\textbf{G}_{j})+\widehat{\textbf{V}}^{*}_{j}\textbf{R}_{\sigma})^{-1}$
$\displaystyle=$
$\displaystyle(\widehat{\textbf{V}}^{*}_{j}\widehat{\textbf{W}}_{j}(\sigma\textbf{I}-\textbf{G}_{j}))^{-1}-(\widehat{\textbf{V}}^{*}_{j}\widehat{\textbf{W}}_{j}(\sigma\textbf{I}-\textbf{G}_{j}))^{-1}\textbf{S}_{j}(\sigma)\widehat{\textbf{V}}_{j}^{*}\textbf{R}_{\sigma}(\widehat{\textbf{V}}^{*}_{j}\widehat{\textbf{W}}_{j}(\sigma\textbf{I}-\textbf{G}_{j}))^{-1}.$

Abbreviating
$\textbf{G}_{\sigma,j}:=\widehat{\textbf{V}}^{*}_{j}\widehat{\textbf{W}}_{j}(\sigma\textbf{I}-\textbf{G}_{j})$
and substituting into (18) yields

$\displaystyle f(\textbf{A})\textbf{b}$ $\displaystyle\approx$
$\displaystyle\widehat{\textbf{V}}_{j}\frac{1}{2\pi
i}\int_{\Gamma}f(\sigma)(\sigma\textbf{I}-\textbf{G}_{j})^{-1}d\sigma\hskip
2.84544pt(\widehat{\textbf{V}}^{*}_{j}\widehat{\textbf{W}}_{j})^{-1}\widehat{\textbf{V}}_{j}^{*}\textbf{b}$
$\displaystyle\mbox{}-\widehat{\textbf{V}}_{j}\frac{1}{2\pi
i}\int_{\Gamma}f(\sigma)\textbf{G}_{\sigma,j}^{-1}(\textbf{I}+\widehat{\textbf{V}}_{j}^{*}\textbf{R}_{\sigma}\textbf{G}_{\sigma,j}^{-1})^{-1}\widehat{\textbf{V}}_{j}^{*}\textbf{R}_{\sigma}\textbf{G}_{\sigma,j}^{-1}\widehat{\textbf{V}}^{*}_{j}\textbf{b}\hskip
2.84544ptd\sigma,$

which is the desired result.

Version 3 of rFOM2 uses numerical quadrature for $\mathcal{I}$ in (19) giving
the approximation

$\displaystyle\widetilde{f}_{3}=$
$\displaystyle\widehat{\textbf{V}}_{j}f(\textbf{G}_{j})(\widehat{\textbf{V}}_{j}^{*}\widehat{\textbf{W}}_{j})^{-1}\widehat{\textbf{V}}_{j}^{*}\textbf{b}$
$\displaystyle-\widehat{\textbf{V}}_{j}\sum_{\ell=1}^{n_{\textrm{quad}}}\omega_{\ell}f(z_{\ell})\textbf{G}_{z_{\ell},j}^{-1}(\textbf{I}+\widehat{\textbf{V}}_{j}^{*}\textbf{R}_{z_{\ell}}\textbf{G}_{z_{\ell},j}^{-1})^{-1}\widehat{\textbf{V}}_{j}^{*}\textbf{R}_{z_{\ell}}\textbf{G}_{z_{\ell},j}^{-1}\widehat{\textbf{V}}_{j}^{*}\textbf{b}.$

The above formulation requires non-singularity of the matrix
$\sigma\textbf{I}-\textbf{G}_{j}$. This is the prime motivation for our
scaling involving D since it allows us to avoid this singularity in the case
when when $\sigma$ is close to $1$. This is not always necessary and in most
cases it should suffice to take $\textbf{D}=\textbf{I}$. We also require the
function $f$ to be defined on the spectrum of D which is not a tight
restriction for most functions of interest. For example, consider choosing D
such that U has unit columns. Then functions with a singularity at zero, would
always be defined on D in practice since we will never have a zero column in
U.

Algorithm 3 rFOM2 \- A recycled FOM approximation for
$f(\textbf{A})\textbf{b}$ (Version 3)

1: Input: $\textbf{A}\in\mathbb{C}^{n\times n}$,
$\textbf{b}\in\mathbb{C}^{n}$, scalar function $f_{s}$, matrix function $f$,
$\textbf{U}\in\mathbb{C}^{n\times k}$ whose columns span the augmenting
recycled subspace $\mathcal{U}$, $\textbf{C}=\textbf{A}\textbf{U}$, Arnoldi
cycle length $j$

2: Build a basis for $\mathcal{K}_{j}(\textbf{A},\textbf{b})$ via the Arnoldi
process, generating $\textbf{V}_{j+1}$ and $\overline{\textbf{H}}_{j}$

3: Choose scaling matrix $\textbf{D}\in\mathbb{C}^{k\times k}$

4: $\textbf{U}\leftarrow\textbf{U}\textbf{D}$,
$\widehat{\textbf{V}}_{j}\leftarrow\begin{bmatrix}\textbf{U}&\textbf{V}_{j}\end{bmatrix}$,
$\widehat{\textbf{W}}_{j}\leftarrow\begin{bmatrix}\textbf{C}&\textbf{V}_{j}\end{bmatrix}$

5: Determine contour containing the spectrum of A and determine quadrature
rule for the integral $\mathcal{I}$ in (19) with nodes $z_{\ell}$ and weights
$\omega_{\ell},\ell=1,\ldots,n_{\textrm{quad}}$

6: Set $\textbf{t}=0\in\mathbb{C}^{j+k}$

7: for $\ell=1,..,n_{\textrm{quad}}$ do

8: $\mu\leftarrow\omega_{\ell}f_{s}(z_{\ell})$

9: $\textbf{t}\leftarrow\textbf{t}+\mu\hskip
2.84544pt\textbf{G}_{z,j}^{-1}(\textbf{I}+\widehat{\textbf{V}}_{j}^{*}\textbf{R}_{z_{\ell}}\textbf{G}_{z,j}^{-1})^{-1}\widehat{\textbf{V}}_{j}^{*}\textbf{R}_{z_{\ell}}\textbf{G}_{z,j}^{-1}\widehat{\textbf{V}}_{j}^{*}\textbf{b}$

10: end for

11:
$\widetilde{f}_{3}\leftarrow\widehat{\textbf{V}}_{j}f(\textbf{G}_{j})(\widehat{\textbf{V}}^{*}_{j}\widehat{\textbf{W}}_{j})^{-1}\widehat{\textbf{V}}^{*}_{j}\textbf{b}-\widehat{\textbf{V}}_{j}\textbf{t}$

12: Update U and C such that $\textbf{C}=\textbf{A}\textbf{U}$ to use for next
matrix A and vector b

An algorithmic description of rFOM2, version 3, is given as Algorithm 3 and an
arithmetic cost summary is given in table Table 3. The cost for
$f(\textbf{G}_{j})$ will depend on whether we can use specifically tailored
methods for specific functions $f$. Hence, we carry a term of
$\mathcal{O}((j+k)^{3})$ for that cost. The cost is proportional to
$(j+k)^{3}$ if instead of a specific method we just use the eigendecomposition
of $\textbf{G}_{j}$. The quantity $\textbf{R}_{z_{\ell}}$ is assumed to be
computed as in version 2.

operation | arithmetic cost | how often
---|---|---
$\widehat{\textbf{V}}_{j}^{*}\widehat{\textbf{W}}_{j}$ | $n(j+k)^{2}$ | 1
Linear solve with $\widehat{\textbf{V}}_{j}^{*}\widehat{\textbf{W}}_{j}$ | $\frac{1}{3}(j+k)^{3}$ | 1
$f(\textbf{G}_{j})$ | $\mathcal{O}((j+k)^{3})$ | 1
$\textbf{U}^{*}(\textbf{U}-\textbf{C})$ | $nk^{2}$ | 1
$\textbf{V}_{j}^{*}(\textbf{U}-\textbf{C})$ | $njk$ | 1
$\widehat{\textbf{V}}_{j}^{*}\widehat{\textbf{W}}_{j}(z\textbf{I}-\textbf{G}_{j}):=\textbf{G}_{z,j}$ | $(j+k)^{3}$ | $n_{quad}$
$\textbf{G}_{z,j}^{-1}$ | $(j+k)^{3}$ | $n_{quad}$
$\widehat{\textbf{V}}_{j}^{*}\textbf{R}_{z_{\ell}}\textbf{G}_{z,j}^{-1}$ | $(j+k)^{3}$ | $n_{quad}$
Linear solve with $(\textbf{I}+\widehat{\textbf{V}}_{j}^{*}\textbf{R}_{z_{\ell}}\textbf{G}^{-1}_{z,j})$ | $\frac{1}{3}(j+k)^{3}$ | $n_{quad}$
total cost:
$n((j+k)^{2}+k^{2}+jk)+n_{\textrm{quad}}\frac{10}{3}(j+k)^{3}+\mathcal{O}((j+k)^{3})$
Table 3: Arithmetic cost of rFOM2 version 3.

We can summarize the cost analysis in this section by saying that, besides the
cost for the Arnoldi process, the costs for all three versions have leading
terms which are $n$ times a sum of quadratic terms in $k$ and $j$ plus
$n_{\textrm{quad}}$ times a sum of cubic terms in $k$ and $j$. In this sense
all versions have comparable cost, although the constants involved in the
quadratic and cubic terms differ from one version to the other. For all
versions, the $n_{quad}$ term in the cost is not dependent on the dominant $n$
and thus one can increase the number of quadrature points without causing
large growths in cost.

##  7 Extension to Stieltjes functions

A Stieltjes function $f$ is a function which admits an integral representation

$f(z)=\int_{-\infty}^{0}\frac{g(\sigma)}{\sigma-z}d\sigma,\enspace
z\in\mathbb{C}\setminus(-\infty,0]$

where $g$ is a function which has constant sign in $[0,\infty).$ Important
Stieltjes functions are the inverse powers $f(z)=z^{-\alpha}$ with
$\alpha\in(0,1)$, as well as other examples such as
$f(z)=\frac{\mbox{log}(1-z)}{z}$: see [13, 22].

The integral representation of a Stieltjes matrix function

$f(\textbf{A})=\int_{-\infty}^{0}g(\sigma)(\sigma\textbf{I}-\textbf{A})^{-1}d\sigma$

resembles the Cauchy integral representation (3) in that it also relies on the
resolvent $(\sigma\textbf{I}-\textbf{A})^{-1}$, the difference being that
instead of a contour $\Gamma$ we now have a fixed integration interval
$(-\infty,0]$. Thus we can also use all three versions of rFOM2 on Stieltjes
functions with the obvious change in the domain of integration. Typical
quadrature rules for the interval $(-\infty,0]$ first map the infinite
interval to a finite interval, $[-1,1)$, say, and then use (variants of)
Gaussian quadrature; see [13], e.g.

##  8 Harmonic Ritz procedure for updating the recycle subspace

A practical implementation of rFOM2 requires an appropriate choice of
augmentation/recycle subspace $\mathcal{U}$. For most functions of practical
interest (eg. log, sign, inverse) the singularity is at the origin. We thus
restrict our attention to these types of functions and aim to recycle the $k$
harmonic Ritz vectors corresponding to the $k$ smallest eigenvalues of the
previous matrix in the sequence.

The augmented Arnoldi relation (17) can be written in its non shifted form as

$\textbf{A}\widehat{\textbf{V}}_{j}=\widehat{\textbf{W}}_{j+1}\overline{\textbf{G}}_{j},$

where
$\widehat{\textbf{W}}_{j+1}=\begin{bmatrix}\textbf{C}&\textbf{V}_{j+1}\end{bmatrix}$
and $\overline{\textbf{G}}_{j}:=\begin{bmatrix}\textbf{D}&\textbf{0}\\\
\textbf{0}&\overline{\textbf{H}}_{j}\end{bmatrix}$.

A procedure for computing harmonic Ritz vectors for the case where the columns
of $\widehat{\textbf{W}}_{j+1}$ are orthonormal has been outlined for the
context of GCRO-DR in [32]. The rFOM2 algorithm is an unprojected augmented
Krylov subspace method and thus we do not have orthonormal columns in
$\widehat{\textbf{W}}_{j+1}$. We present a slight modification of the harmonic
Ritz procedure outlined in [32] for this case.

######  Proposition 8.1.

The harmonic Ritz problem for the matrix A with respect to the augmented
Krylov subspace space $\mathcal{K}_{j}(\textbf{A},\textbf{b})+\mathcal{U}$
involves solving the following eigenproblem

$\overline{\textbf{G}}^{*}_{j}\widehat{\textbf{W}}^{*}_{j+1}\widehat{\textbf{W}}_{j+1}\overline{\textbf{G}}_{j}\textbf{g}_{i}=\theta_{i}\overline{\textbf{G}}_{j}^{*}\widehat{\textbf{W}}^{*}_{j+1}\widehat{\textbf{V}}_{j}\textbf{g}_{i}.$

######  Proof 8.2.

The harmonic Ritz problem for the matrix A with respect to the space
$\mathcal{K}_{j}(\textbf{A},\textbf{b})+\mathcal{U}$ requires one to solve the
following problem

$\displaystyle\mbox{Find}\hskip 14.22636pt(\textbf{y}_{i},\mu_{i})\hskip
14.22636pt\mbox{such that}\hskip
14.22636pt\textbf{A}^{-1}\textbf{y}_{i}-\mu_{i}\textbf{y}_{i}\perp\textbf{A}(\mathcal{K}_{j}(\textbf{A},\textbf{b})+\mathcal{U})$
$\displaystyle\hskip 14.22636pt\mbox{with}\hskip
14.22636pt\textbf{y}_{i}\in\textbf{A}(\mathcal{K}_{j}(\textbf{A},\textbf{b})+\mathcal{U})$

and set $\theta_{i}=1/\mu_{i}$.

Since $\textbf{y}_{i}$ is taken from the space
$\textbf{A}(\mathcal{K}_{j}(\textbf{A},\textbf{b})+\mathcal{U})$, it can be
written as $\textbf{y}_{i}=\textbf{A}\widehat{\textbf{V}}_{j}\textbf{g}_{i}$
for some $\textbf{g}_{i}\in\mathbb{C}^{j+k}$. Application of the orthogonality
constraint yields an eigenproblem for $\textbf{g}_{i}$

$\displaystyle
0=(\textbf{A}\widehat{\textbf{V}}_{j})^{*}(\textbf{A}^{-1}\textbf{y}_{i}-\mu_{i}\textbf{y}_{i})=(\textbf{A}\widehat{\textbf{V}}_{j})^{*}(\textbf{A}^{-1}\textbf{A}\widehat{\textbf{V}}_{j}\textbf{g}_{i}-\mu_{i}\textbf{A}\widehat{\textbf{V}}_{j}\textbf{g}_{i})$
$\displaystyle=(\widehat{\textbf{W}}_{j+1}\overline{\textbf{G}}_{j})^{*}(\widehat{\textbf{V}}_{j}\textbf{g}_{i}-\mu_{i}\widehat{\textbf{W}}_{j+1}\overline{\textbf{G}}_{j}\textbf{g}_{i})=\overline{\textbf{G}}^{*}_{j}\widehat{\textbf{W}}^{*}_{j+1}\widehat{\textbf{V}}_{j}\textbf{g}_{i}-\mu_{i}\overline{\textbf{G}}_{j}^{*}\widehat{\textbf{W}}^{*}_{j+1}\widehat{\textbf{W}}_{j+1}\overline{\textbf{G}}_{j}\textbf{g}_{i}$
$\displaystyle\implies\overline{\textbf{G}}^{*}_{j}\widehat{\textbf{W}}^{*}_{j+1}\widehat{\textbf{W}}_{j+1}\overline{\textbf{G}}_{j}\textbf{g}_{i}=\theta_{i}\overline{\textbf{G}}_{j}^{*}\widehat{\textbf{W}}^{*}_{j+1}\widehat{\textbf{V}}_{j}\textbf{g}_{i}.$

We then set $\textbf{y}_{i}=\textbf{A}\widehat{\textbf{V}}_{j}\textbf{g}_{i}$.

The difference between this problem and that of the harmonic Ritz problem in
[32] is the inclusion of the factor
$\widehat{\textbf{W}}^{*}_{j+1}\widehat{\textbf{W}}_{j+1}$. We also note since
rFOM2 does not require orthonormal columns of C, there is no need for an
additional QR factorization after computing $\textbf{C}=\textbf{A}\textbf{U}$.

##  9 Numerical experiments

In this section we present results of our numerical experiments. These
experiments aim to test the effectiveness of all three implementations of
rFOM2 as both an augmented Krylov subspace method, and as a method for
recycling. Our experiments require knowledge of an exact solution for a
sequence of slowly changing matrices, and have thus been limited to relatively
small problem sizes. We also note that the methods are sensitive to loss of
orthogonality in the Arnoldi vectors, and thus propose to reorthogonalize
these vectors to maintain orthogonality to within a good level of accuracy.

We compare the results of rFOM2 to the standard Arnoldi approximation (6). For
the evaluation of the matrix function on the small $j\times j$ and
$(j+k)\times(j+k)$ matrices in standard Arnoldi and version 3 (line 11 in
Algorithm 3), respectively, we use a “direct” method, given by an appropriate
built-in Matlab function. These MATLAB functions also allow us to obtain what
we considered the “exact” solution for each experiment which we used to
determine the error. Our experiments, in addition, compute the Arnoldi
approximation evaluated using the same quadrature procedure on the Hessenberg
matrix as rFOM2 (denoted Arnoldi (q)).

All experiments involving the sign function have been performed using a
$4^{4}$ lattice QCD matrix of size $3072\times 3072$ given by
$\textbf{Q}=\Gamma_{5}\textbf{D}$, where D is the Wilson Dirac matrix
available as matrix ID $1591$ from the SuiteSparse Matrix Collection [8]. The
sign function is typically computed using the relation
$\mbox{sign}(\textbf{Q})=(\textbf{Q}^{2})^{-1/2}\textbf{Q}$, and our
experiments actually report the results for the inverse square root of
$\textbf{Q}^{2}$, which can be represented as a Stieltjes integral (see
Section 7)

$z^{-1/2}=\frac{-1}{\pi}\int_{-\infty}^{0}\frac{\sigma^{-1/2}}{\sigma-z}d\sigma,$

for which we use a variant of Gaussian quadrature as detailed in [13].

For all other functions tested we take the contour of integration $\Gamma$ to
be a circle of radius $r$ centered at the point $c$ and use the trapezoidal
rule for numerical integration. The contour $\Gamma$ is given by the
parameterization

$z=re^{i\theta}+c,\hskip 8.5359pt0\leq\theta\leq 2\pi,$

and the nodes of the trapezoidal rule are the $n_{\textrm{quad}}$ equidistant
points given by

$z_{\ell}=re^{i\theta_{\ell}}+c\hskip 2.84544pt,\hskip
2.84544pt\theta_{\ell}=\ell\Delta\theta,\hskip
2.84544pt\ell=0,1,..,n_{quad}-1,\enspace\Delta\theta=2\pi/n_{quad}.$

Since $\Gamma$ is closed, the trapezoidal rule is particularly accurate and
features exponential convergence; see [39, 41], e.g.

###  9.1 Results

(a)

(b)

(c)

(d)

Figure 1:  Comparison of Arnoldi and the three versions of rFOM2 for computing
$f(\textbf{A})\textbf{b}.$ QCD examples: (a) $f(z)=z^{-\frac{1}{2}}$ and
$\textbf{A}=\textbf{Q}^{2}$, (b) $f(z)=\frac{1}{z}$ and
$\textbf{A}=\textbf{D}$. Finite difference examples: (c) $f(z)=\mbox{log}(z)$
and A a 2d convection diffusion matrix (d) $f(z)=\sqrt{z}$ and A discrete
Laplace matrix. All experiments use $j=40$ Arnoldi iterations and recycle
space dimension $k=20$.

Figure 1 shows the results obtained from augmenting the Krylov subspace with
high quality eigenvector approximations computed directly from the matrix A.
We compared the three versions of rFOM2 (denoted
$\widetilde{f}_{1},\widetilde{f}_{2}$ and $\widetilde{f}_{3}$ respectively) to
the the standard Arnoldi approximation with “direct” evaluation of $f$ on the
Hessenberg matrix (Arnoldi) and the quadrature based approximation (Arnoldi
(q)). The experiments were run for a range of values for $n_{quad}$.

The experiment in Figure 1 (a) corresponds to the inverse square root of
$\textbf{Q}^{2}$. In Figure 1 (b) we tested the method on the inverse function
of the standard Wilson Dirac matrix D. Both these experiments show that all
implementations of rFOM2 lead to a significant reduction in the error as
compared to the Arnoldi approximation. The experiments in Figure 1 (c/d)
further demonstrate the effectiveness of rFOM2 using rather academic examples
for which the standard Arnoldi approximation performs poorly, while rFOM2
significantly increases the quality of the approximation. In (c) we compute
the log function of a matrix arising from a second order finite difference
discretization of a 2D convection diffusion problem on a $100\times 100$ mesh,
taken from the open source funm$\\_$quad code [13]. In (d) we compute the
square root of the discrete Laplace matrix corresponding to a discretization
of Poisson’s equation on the unit square as given in the MATLAB test matrix
gallery.

There are some common observations we can make from these experiments.
Firstly, we observe that beyond a certain number of quadrature points, the
error of quadrature Arnoldi stagnates and becomes identical to the error of
standard Arnoldi. This indicates that the quadrature error has become
negligible as compared to the subspace approximation error. For the augmented
methods, one typically needs more quadrature points before such a stagnation
occurs.

At stagnation and beyond, the rFOM2 approximations significantly outperform
the standard Arnoldi approximation and all three implementations produce the
same approximation. We also note that the $\widetilde{f}_{3}$ approximation is
more accurate than the other two versions of rFOM2 for the same number of
quadrature points. This is because this approximation involves a term which
does not depend on quadrature but a “closed form” evaluation of the function.
The other two approximations are based entirely on quadrature and thus require
a larger number of quadrature points to achieve the same level of accuracy.
Even though time measurements cannot be regarded as conclusive with our
implementations and software, we just mention that versions $2$ and $3$ tend
to take somewhat longer compute time than version $1$, an observation which is
in line with our cost analysis.

In Figure 2, 3 and 4 we demonstrate the effectiveness of the three versions of
rFOM2 as a recycle method. For each experiment we use enough quadrature points
so that all three implementations produce the same approximation. We denote
the common error curve for all implementations as rFOM2 in the plots.

(a)

(b)

Figure 2: Comparison of standard Arnoldi and the three versions of rFOM2 when
applying the inverse of the $4^{4}$ lattice QCD matrix D on a sequence of 10
random vectors. All experiments used $j=50$ Arnoldi iterations, a recycle
subspace dimension $k=20$ and $1000$ quadrature points. $\epsilon=0$ in (a),
$\epsilon=10^{-3}$ in (b).

The experiments in Figure 2 involve computing the inverse function on a
sequence of slowly changing Wilson Dirac matrices against 10 random vectors.
We used a parameter $\epsilon$ to describe the strength of the perturbation
between each matrix in the sequence. In Figure 2 (a) we demonstrate the result
for $\epsilon=0$ which corresponds to the special case where the matrix
$\textbf{A}^{(i)}$ remains fixed for all $i$, and only the vectors
$\textbf{b}^{(i)}$ change throughout the sequence. In Figure 2 (b) we
demonstrate the case where the matrices now change according to a random
perturbation of strength $\epsilon=10^{-3}$. For all our experiments we make
no assumption on the relationship between consecutive vectors
$\textbf{b}^{(i)}$ and generate an independent random vector for each new
problem in the sequence. The subspace $\mathcal{U}$ has been computed between
each system using the harmonic Ritz procedure outlined in Section 8. In the
plots below, enough quadrature points have been chosen so that the standard
Arnoldi and the quadrature based Arnoldi produce the same approximation, and
are thus represented as the same curve. From these experiments we see that all
three variants of rFOM2 exhibit significant reductions in error as compared to
the standard Arnoldi approximation without recycling, and that they
consistently reduce the error further as we move throughout the sequence of
problems. This is one of the key qualities of an effective recycling method as
outlined in [32].

(a)

(b)

(c)

(d)

Figure 3: Comparison of Arnoldi approximation with the 3 versions of rFOM2 for
a sequence of 30 applications of the inverse square root function with
randomly changing matrices. The strength of change is controlled by
$\epsilon$. Fig (a) is for $\epsilon=0$ (matrices remain constant), and
$\epsilon=10^{-4}$ in (b), $\epsilon=10^{-3}$ in (c), and $\epsilon=10^{-2}$
in (d). All experiments use $j=50$ Arnoldi iterations, recycle subspace
dimension of size $k=20$ and $n_{\textrm{quad}}=30$.

In Figure 3 we repeat a similar experiment as in Figure 2, but this time with
30 computations of $f(\textbf{A}^{(i)})\textbf{b}^{(i)}$ with $f$ the inverse
square root, $\textbf{A}^{(1)}=\textbf{Q}^{2}$ and then producing a sequence
of perturbed Hermitian matrices $\textbf{A}^{(i)}$ with the strength of change
controlled by $\epsilon$ as in the previous experiment. Figure 3 (a/b/c) were
performed with $\epsilon=0$, $10^{-4}$ and $10^{-3}$, respectively, and
demonstrate the same behaviour as observed in the previous experiment. The
error is reduced considerably as we move along the sequence of matrix function
applications.

Figure 4: Minimum eigenvalue for each $\textbf{A}^{(i)}$ in the experiments in
Figure 3 for the different values of $\epsilon$. For $\epsilon=0.01$, negative
eigenvalues appear from the $14^{th}$ system on.

The plot in Figure 3 (d) has $\epsilon=10^{-2}$. This experiment performs
similarly to the other experiments up until problem index 14, where the random
perturbations have caused the spectrum to shift to the negative real part of
the complex plane (see Figure 4). This highlights an important practical
limitation of rFOM2: The contour of integration (which we consider to be the
interval $(-\infty,0]$ in the Stieltjes case) needs to enclose the spectrum of
the matrix. If this is not the case anymore, as in problem number 14, where A
has a negative eigenvalue, we do not compute the correct integral unless we
adapt the contour $\Gamma$. This issue will be highly application dependent.
It will not pose an issue for applications where sets containing the spectra
of all matrices involved are known in advance or for the special case where
the matrix remains fixed.

Figure 5: Accuracy of the recycle subspace $\mathcal{U}$ as an eigenvector
approximation from the experiments in Figure 3.

Finally, we demonstrate how the recycle subspace $\mathcal{U}$ converges to a
good approximate eigenspace for each of the problems in the sequence. This is
one of the requirements for an effective recycled method outlined in [32]. We
demonstrate this in Figure 5 using the recycle subspace $\mathcal{U}$ for each
of the experiments in Figure 3. Figure 5 shows how accurately the subspace
$\mathcal{U}$ approximates the eigenvectors of A as the sequence of matrix
function applications progresses. This was measured by computing the angle
$\theta(\mathcal{U},\mathcal{Z})$ between the subspace $\mathcal{U}$ and the
subspace $\mathcal{Z}$ spanned by the $k$ smallest exact eigenvectors of each
matrix $\textbf{A}^{(i)}$ in the sequence. In general, we see that the quality
of the computed recycled subspace improves, i.e.
$\theta(\mathcal{U},\mathcal{Z})$ becomes smaller, as the sequence progresses
provided the matrices in the sequence are slowly changing. For a large enough
perturbation of the matrix, the quality of $\mathcal{U}$ appears to remain
roughly fixed throughout the sequence.

##  10 Conclusions

In this paper we have presented a new algorithm rFOM2 for computing an
augmented Krylov subspace approximation to the action of a matrix function on
a vector. The approximation allows one to recycle information between
sequences of matrix function applications wherein the matrix and/or right hand
side is changing. We presented three different approaches to an implementation
and have demonstrated their effectiveness using numerical experiments. Our
experiments suggest the method works very well as a recycling method, but has
potential limitations when sufficient a priori information of each matrix in
the sequence is not available This is not a major issue for most applications
of interest such as Lattice QCD where we know where the spectrum is contained.

## Acknowledgments

The first author would like to thank The Hamilton Scholars for funding this
work, and Kathryn Lund for useful conversations.

## References

  * [1] G. Arnold, N. Cundy, J. v. d. Eshof, A. Frommer, S. Krieg, T. Lippert, and K. Schäfer, Numerical methods for the QCD overlap operator: Ii. Optimal Krylov subspace methods, in QCD and Numerical Analysis III, Springer, 2005, pp. 153–167.
  * [2] J. Bloch, A. Frommer, B. Lang, and T. Wettig, An iterative method to compute the sign function of a non-Hermitian matrix and its application to the overlap Dirac operator at nonzero chemical potential, Computer Physics Communications, 177 (2007), pp. 933–943.
  * [3] A. Boriçi, A. Frommer, B. Joó, A. Kennedy, and B. Pendleton, QCD and Numerical Analysis iii, in Proceedings of the third international workshop on Numerical Analysis and Lattice QCD, Edinburgh, UK, Springer, 2003.
  * [4] J. Brannick, A. Frommer, K. Kahl, B. Leder, M. Rottmann, and A. Strebel, Multigrid preconditioning for the overlap operator in Lattice QCD, Numer. Math., 132 (2016), pp. 463–490.
  * [5] L. Burke and K. M. Soodhalter, Augmented unprojected Krylov subspace methods from an alternative view of an existing framework, arXiv preprint arXiv:2206.12315, (2022).
  * [6] N. Cundy, S. Krieg, G. Arnold, A. Frommer, T. Lippert, and K. Schilling, Numerical methods for the QCD overlap operator iv: Hybrid Monte Carlo, Computer Physics Communications, 180 (2009), pp. 26–54.
  * [7] N. Cundy, J. van den Eshof, A. Frommer, S. Krieg, T. Lippert, and K. Schäfer, Numerical methods for the QCD overlap operator: Iii. Nested iterations, Computer Physics Communications, 165 (2005), pp. 221–242.
  * [8] T. A. Davis and Y. Hu, The University of Florida sparse matrix collection, ACM TOMS, 38 (2011), pp. 1–25, https://doi.org/10.1145/2049662.2049663.  Collection available through https://sparse.tamu.edu/about.
  * [9] E. de Sturler, Nested Krylov methods based on GCR, Journal of Computational and Applied Mathematics, 67 (1996), pp. 15–41.
  * [10] M. Eiermann and O. G. Ernst, A restarted Krylov subspace method for the evaluation of matrix functions, SIAM Journal on Numerical Analysis, 44 (2006), pp. 2481–2504.
  * [11] M. Eiermann, O. G. Ernst, and S. Güttel, Deflated restarting for matrix functions, SIAM Journal on Matrix Analysis and Applications, 32 (2011), pp. 621–641.
  * [12] A. Frommer, S. Güttel, and M. Schweitzer, Convergence of restarted Krylov subspace methods for Stieltjes functions of matrices, SIAM Journal on Matrix Analysis and Applications, 35 (2014), pp. 1602–1624.
  * [13] A. Frommer, S. Güttel, and M. Schweitzer, Efficient and stable Arnoldi restarts for matrix functions based on quadrature, SIAM Journal on Matrix Analysis and Applications, 35 (2014), pp. 661–683, https://doi.org/10.1137/13093491X.
  * [14] A. Frommer, T. Lippert, B. Medeke, and K. Schilling, Numerical Challenges in Lattice Quantum Chromodynamics: Joint Interdisciplinary Workshop of John Von Neumann Institute for Computing, Jülich, and Institute of Applied Computer Science, Wuppertal University, August 1999, vol. 15, Springer Science & Business Media, 2000.
  * [15] A. Frommer, K. Lund, and D. B. Szyld, Block Krylov subspace methods for functions of matrices, Electronic Transactions on Numerical Analysis, 47 (2017), pp. 100–126.
  * [16] A. Frommer, K. Lund, and D. B. Szyld, Block Krylov subspace methods for functions of matrices ii: Modified block FOM, SIAM Journal on Matrix Analysis and Applications, 41 (2020), pp. 804–837.
  * [17] A. Gaul, Recycling Krylov subspace methods for sequences of linear systems, PhD thesis, Technische Universität Berlin, Fakultät II - Mathematik und Naturwissenschaften, 2014.
  * [18] A. Gaul, M. H. Gutknecht, J. Liesen, and R. Nabben, A framework for deflated and augmented Krylov subspace methods, SIAM Journal on Matrix Analysis and Applications, 34 (2013), pp. 495–518.
  * [19] M. H. Gutknecht, Spectral deflation in Krylov solvers: A theory of coordinate space based methods, Electronic Transactions on Numerical Analysis, 39 (2012), pp. 156–185.
  * [20] M. H. Gutknecht, Deflated and augmented Krylov subspace methods: A framework for deflated BiCG and related solvers, SIAM Journal on Matrix Analysis and Applications, 35 (2014), pp. 1444–1466.
  * [21] N. Hale, N. J. Higham, and L. N. Trefethen, Computing $a^{\alpha}$, $\mathrm{log}(a)$, and related matrix functions by contour integrals, SIAM Journal on Numerical Analysis, 46 (2008), pp. 2505–2523.
  * [22] P. Henrici, Applied and computational complex analysis. Vol. 2, Wiley Classics Library, John Wiley & Sons, Inc., New York, 1991.
  * [23] N. J. Higham, Functions of Matrices, SIAM, Philadelphia, 2008.
  * [24] M. Hochbruck and C. Lubich, On Krylov subspace approximations to the matrix exponential operator, SIAM Journal on Numerical Analysis, 34 (1997), pp. 1911–1925.
  * [25] M. Hochbruck and A. Ostermann, Exponential Integrators, Acta Numerica, 19 (2010), pp. 209–286.
  * [26] R. A. Horn and C. R. Johnson, Topics in matrix analysis, Cambridge University Press, Cambridge, 1994.  Corrected reprint of the 1991 original.
  * [27] B. Joó, C. Jung, N. H. Christ, W. Detmold, R. G. Edwards, M. Savage, and P. Shanahan, Status and future perspectives for lattice gauge theory calculations to the exascale and beyond, The European Physical Journal A, 55 (2019), pp. 1–26.
  * [28] F. Knechtli, M. Günther, and M. Peardon, Lattice Quantum Chromodynamics: Practical Essentials, Springer, 2017.
  * [29] K. Lund, A new block Krylov subspace framework with applications to functions of matrices acting on multiple vectors, PhD thesis, Department of Mathematics, Temple University, 2018.
  * [30] R. B. Morgan, GMRES with deflated restarting, SIAM Journal on Scientific Computing, 24 (2002), pp. 20–37.
  * [31] H. Neuberger, Exactly massless quarks on the lattice, Physics Letters B, 417 (1998), pp. 141–144.
  * [32] M. L. Parks, E. De Sturler, G. Mackey, D. D. Johnson, and S. Maiti, Recycling Krylov subspaces for sequences of linear systems, SIAM Journal on Scientific Computing, 28 (2006), pp. 1651–1674.
  * [33] Y. Saad, Iterative methods for sparse linear systems, SIAM, 2003.
  * [34] V. Simoncini, Restarted full orthogonalization method for shifted linear systems, BIT Numerical Mathematics, 43 (2003), pp. 459–466.
  * [35] B. Singer and S. Spilerman, The representation of social processes by Markov models, American journal of sociology, 82 (1976), pp. 1–54.
  * [36] K. M. Soodhalter, Block Krylov subspace recycling for shifted systems with unrelated right-hand sides, SIAM Journal on Scientific Computing, 38 (2016), pp. A302–A324.
  * [37] K. M. Soodhalter, E. de Sturler, and M. E. Kilmer, A survey of subspace recycling iterative methods, GAMM‐Mitteilungen, 43 (2020).
  * [38] K. M. Soodhalter, D. B. Szyld, and F. Xue, Krylov subspace recycling for sequences of shifted linear systems, Applied Numerical Mathematics, 81 (2014), pp. 105–118.
  * [39] L. N. Trefethen and J. A. C. Weideman, The exponentially convergent trapezoidal rule, SIAM Review, 56 (2014), pp. 385–458.
  * [40] J. van den Eshof, A. Frommer, T. Lippert, K. Schilling, and H. A. van der Vorst, Numerical methods for the QCDd overlap operator. i. Sign-function and error bounds, Computer Physics Communications, 146 (2002), pp. 203–224.
  * [41] J. A. C. Weideman, Numerical integration of periodic functions: A few examples, The American Mathematical Monthly, 109 (2002), pp. 21–36.