in feature removal.

###  F.2 Empirical observations

In eq. F.1.1, the sampling-dependent term
$g_{\text{sampling}}(\text{removal})$ does not depend on the input
perturbation strength $\lVert x-x^{\prime}\rVert_{2}$, so we empirically
estimate $g_{\text{sampling}}(\text{removal})$ with the intercept of the
empirical upper bound (as illustrated in Figure 3) when $\lVert
x-x^{\prime}\rVert_{2}=0$. As shown in Figure 13 (left),
$g_{\text{sampling}}(\textbf{baseline})$ is indeed zero, whereas
$g_{\text{sampling}}(\textbf{marginal})$ and
$g_{\text{sampling}}(\textbf{conditional})$ empirically decrease with
increasing sample size for feature removal as expected. We observe that the
sampling error tends to be lower for the conditional approach compared to the
marginal approach, potentially due to lower variances in the conditional
versus marginal distributions. Finally, Figure 13 (right) confirms that the
perturbation-dependent term $g(\text{removal})\cdot
h(\text{summary})\cdot\lVert x-x^{\prime}\rVert_{2}$ does not depend on the
removal sample size.

Figure 13: Intercepts and slopes of empirical and theoretical upper bounds
with respect to the input perturbation norm, for Occlusion, Banzhaf, Shapley,
RISE, and LIME attribution differences in the synthetic experiment in Section
A.1, and with varying sample size for estimating feature removal.

##  Appendix G Feature grouping

Rather than generate attribution scores for each feature $x_{i}$, another
option is to partition the feature set and generate group-level scores. This
is common in practice for images [53, 12, 16], where removing individual
pixels is unlikely to produce meaningful prediction differences, and reducing
the number of features makes attributions less computationally intensive. For
example, LIME pre-computes superpixels using a basic segmentation algorithm
[53], and ViT Shapley uses grid-shaped patches defined by a vision transformer
[16].

When computing group-level attributions, we require a set of feature index
groups that we denote as $\\{G_{1},\ldots,G_{g}\\}$, where each group
satisfies $G_{i}\subseteq[d]$. We assume that the $g$ groups are non-
overlapping and cover all the indices, so that we have $G_{i}\cap
G_{j}=\\{\\}$ for $i\neq j$ and $G_{1}\cup\ldots\cup G_{g}=[d]$. Under this
approach, the attributions are a vector of the form
$\phi(f,x)\in\mathbb{R}^{g}$, and we compute them by querying the model only
with subsets that are unions of the groups. For this, we let $T\subseteq[g]$
denote a subset of group indices and denote the corresponding union of groups
as $G_{T}=\bigcup_{i\in T}G_{i}\subseteq[d]$. Using this notation and taking
the Shapley value as an example, the attribution for the $i$th group $G_{i}$
is defined as

$\phi_{i}(f,x)=\frac{1}{g}\sum_{T\subseteq[g]\setminus\\{i\\}}\binom{g-1}{|T|}^{-1}\left(f(x_{G_{T\cup\\{i\\}}})-f(x_{G_{T}})\right),$

where $f(x_{G_{T}})$ is defined based on the chosen feature removal approach.

To generalize our earlier results to this setting, we first remark that
querying the model with each subset induces a vector $v\in\mathbb{R}^{2^{g}}$
containing the predictions for each union $G_{T}$ with $T\subseteq[g]$. The
attributions are then calculated according to $\phi(f,x)=Av$, where
$A\in\mathbb{R}^{g\times 2^{g}}$ is the linear operator associated with each
summary technique (see Table 3). Our approach for bounding the attribution
difference under either input or model perturbation is similar here, because
we can use the same inequalities from Lemma 7, or
$||A(v-v^{\prime})||_{2}\leq||A||_{2}\cdot||v-v^{\prime}||_{2}$ and
$||A(v-v^{\prime})||_{2}\leq||A||_{1,\infty}\cdot||v-v^{\prime}||_{\infty}$,
depending on which norm is available for the vectors of predictions
$||v-v^{\prime}||$.

Following the same proof techniques used previously, we can guarantee the
following bounds in the feature grouping setting.

Input perturbation. Lemma 1 and Lemma 3 guarantee that $f(x_{S})$ is Lipschitz
continuous for any subset $S\subseteq[d]$. These results automatically apply
to the specific subsets induced by feature grouping. If $f(x_{S})$ is
$L^{\prime}$-Lipschitz continuous for all $S\subseteq[d]$, we therefore have
the following inequality,

$|f(x_{G_{T}})-f(x^{\prime}_{G_{T}})|\leq
L^{\prime}\cdot||x_{G_{T}}-x^{\prime}_{G_{T}}||_{2},$

for all group subsets $T\subseteq[g]$. As a result, we can guarantee that the
induced vectors $v,v^{\prime}\in\mathbb{R}^{2^{g}}$ for two inputs
$x,x^{\prime}$ have the following bound:

$||v-v^{\prime}||_{\infty}\leq L^{\prime}\cdot||x-x^{\prime}||_{2}.$

Finally, we see that Theorem 1 only needs to be changed by replacing $d$ with
$g$ in the $h(\text{summary})$ term. This bound may even be tighter when we
use feature groups, because the Lipschitz constant may be smaller when we
restrict our attention to a smaller number of subsets.

Model perturbation. Lemma 5 and Lemma 6 guarantee that the predictions
$f(x_{S})$ and $f^{\prime}(x_{S})$ for two models $f,f^{\prime}$ can differ by
no more than the functional difference $||f-f^{\prime}||$, where the specific
norm depends on the feature removal approach. These results apply to all
subsets $S\subseteq[d]$, so they automatically apply to those induced by
feature grouping. We therefore have the following inequality,

$|f(x_{G_{T}})-f^{\prime}(x_{G_{T}})|\leq||f-f^{\prime}||,$

for all $x$ and all group subsets $T\subseteq[g]$. For the induced vectors
$v,v^{\prime}\in\mathbb{R}^{2^{g}}$, we also have

$||v-v^{\prime}||_{\infty}\leq||f-f^{\prime}||.$

Finally, we see that Theorem 2 applies as well, where the only change required
is replacing $d$ with $g$ in the $h(\text{summary})$ term.

##  Appendix H LIME weighted least squares linear operator

In this section, we discuss the linear operator $A\in\mathbb{R}^{d\times
2^{d}}$ associated with LIME [53]. We begin with a theoretical
characterization that focuses on LIME’s default implementation choices, and we
then corroborate these points with empirical results.

###  H.1 Theory

The observation that LIME is equivalent to a linear operator
$A\in\mathbb{R}^{d\times 2^{d}}$ is discussed in Appendix D. We pointed out
there that LIME’s specific operator depends on several implementation choices,
including the choice of weighting kernel, regularization and intercept term.
Here, we consider that no regularization is used, that an intercept term is
fit, and we focus on the default weighting kernel in LIME’s popular
implementation.333https://github.com/marcotcr/lime The default is the
following exponential kernel,

$w(S)=\exp\left(-\frac{\left(1-\sqrt{|S|/d}\right)^{2}}{\sigma^{2}}\right),$
(H.1.1)

where $\sigma^{2}>0$ is a hyperparameter. The user can set this parameter
arbitrarily, but LIME also considers default values for each data type:
$\sigma^{2}=0.25$ for images, $\sigma^{2}=25$ for language, and
$\sigma^{2}=0.75\sqrt{d}$ for tabular data. We find that these choices
approach two limiting cases, $\sigma^{2}\to\infty$ and $\sigma^{2}\to 0$,
where LIME reduces to other summary techniques; namely, we find that these
limiting cases respectively recover the Banzhaf value [9] and leave-one-out
[67] summary techniques. We argue this first from a theoretical perspective,
and Section H.2 then provides further empirical evidence.

Limiting case $\bm{\sigma^{2}\to\infty}$. In this case, it is easy to see from
eq. H.1.1 that we have $w(S)\to 1$ for all $S\subseteq[d]$. The weighted least
squares objective for LIME then reduces to an unweighted least squares
problem, which according to results from Hammer and Holzman [32] recovers the
Banzhaf value. We might expect that the default settings of $\sigma^{2}$
approach this outcome in practice (e.g., $\sigma^{2}=25$), and we verify this
in our empirical results.

Limiting case $\bm{\sigma^{2}\to 0}$. As we let $\sigma^{2}$ approach $0$,
increasing weight is applied to $w(S)$ when the cardinality $|S|$ is large.
Similar to a softmax activation, the differences between cardinalities become
exaggerated with low temperatures; for very low temperatures, the weight
applied to $|S|=d$ is much larger than the weight applied to $|S|=d-1$, which
is much larger than the weight for $|S|=d-2$, etc. In this regime, applying
weight only to the subset $S=[d]$ leaves the problem undetermined, so the
small residual weight on $|S|=d-1$ plays an important role. The problem begins
to resemble one where only subsets with $|S|\geq d-1$ are taken into account,
which yields a unique solution: following Covert et al. [17], this reduces
LIME to leave-one-out (Occlusion). More formally, if we denote
$A_{\text{LIME}}(\sigma^{2})$ as the linear operator resulting from a specific
$\sigma^{2}$ value, we argue that $\lim_{\sigma^{2}\to
0}A_{\text{LIME}}(\sigma^{2})$ is the operator associated with leave-one-out.
We do not prove this result due to the difficulty of deriving
$A_{\text{LIME}}(\sigma^{2})$ for a specific $\sigma^{2}$ value, but Section
H.2 provides empirical evidence for this claim.

###  H.2 Empirical analysis

Here, we denote the LIME linear operator with the exponential kernel as
$A_{\text{LIME}}$ instead of $A_{\text{LIME}}(\sigma^{2})$ for brevity. Recall
from Appendix D that the LIME linear operator has the form
$A_{\text{LIME}}=(Z^{\prime\top}WZ^{\prime})^{-1}_{[1:]}Z^{\prime\top}W$,
where $Z^{\prime}=\\{0,1\\}^{2^{d}\times(d+1)}$ is a binary matrix with all
ones in the first column and the other columns indicating whether each feature
is in each subset, and $W\in{\mathbb{R}}^{2^{d}\times 2^{d}}$ is a diagonal
matrix with $W_{S,S}=w(S)$ for all $S\in[d]$ and zeros elsewhere.

To empirically compare the linear operator of LIME to those of the Banzhaf
value and leave-one-out (Occlusion), we first derive closed-form expressions
of $Z^{\prime\top}WZ^{\prime}$ and $Z^{\prime\top}W$ as in the following
proposition.

######  Proposition 2.

Define

$\displaystyle a_{0}$
$\displaystyle\equiv\sum_{k=0}^{d}\binom{d}{k}\exp\left(-\left(1-\sqrt{k/d}\right)^{2}/\sigma^{2}\right),$
$\displaystyle a_{1}$
$\displaystyle\equiv\sum_{k=1}^{d}\binom{d-1}{k-1}\exp\left(-\left(1-\sqrt{k/d}\right)^{2}/\sigma^{2}\right),\text{
and}$ $\displaystyle a_{2}$
$\displaystyle\equiv\sum_{k=2}^{d}\binom{d-2}{k-2}\exp\left(-\left(1-\sqrt{k/d}\right)^{2}/\sigma^{2}\right).$

Then $(Z^{\prime\top}WZ^{\prime})$ is a $(d+1)\times(d+1)$ matrix with the
following form:

$\displaystyle(Z^{\prime\top}WZ^{\prime})=\begin{pmatrix}a_{0}&B\\\
B^{\top}&C\end{pmatrix},$

where $B=(a_{1},...,a_{1})\in{\mathbb{R}}^{1\times d}$ , and $C$ is a $d\times
d$ matrix with $a_{1}$ on the diagonal entries and $a_{2}$ on the off-diagonal
entries. Furthermore, for all $S\in[d]$

$(Z^{\prime\top}W)_{iS}=\begin{cases}w(S),&\text{if $i=1$}\\\
\mathbbm{1}\\{i-1\in S\\}w(S)&\text{otherwise}.\end{cases}$

###### Proof.

We first find an expression for each element in $Z^{\prime\top}WZ^{\prime}$.
Consider the first diagonal entry:

$\displaystyle(Z^{\prime\top}WZ^{\prime})_{11}=Z^{\prime\top}_{1}WZ^{\prime}_{1}$
$\displaystyle={\bm{1}}_{2^{d}}^{\top}W{\bm{1}}_{2^{d}}$
$\displaystyle=\sum_{S\subseteq[d]}w(S)$
$\displaystyle=\sum_{k=0}^{d}\binom{d}{k}\exp\left(-\left(1-\sqrt{k/d}\right)^{2}/\sigma^{2}\right)=a_{0}.$

For the rest of the diagonal entries, we have:

$\displaystyle(Z^{\prime\top}WZ^{\prime})_{ii}=Z^{\prime\top}_{i}WZ^{\prime}_{i}$
$\displaystyle=\sum_{S:i-1\in S}w(S)$
$\displaystyle=\sum_{k=1}^{d}\binom{d-1}{k-1}\exp\left(-\left(1-\sqrt{k/d}\right)^{2}/\sigma^{2}\right)=a_{1},$

for $i=2,...,d+1$. For the off-diagonal entries, we have:

$(Z^{\prime\top}WZ^{\prime})_{1i}=(Z^{\prime\top}WZ^{\prime})_{i1}={\bm{1}}_{2^{d}}^{\top}WZ^{\prime}_{i}=\sum_{S:i-1\in
S}w(S)=a_{1},$

for $i=2,...,d+1$. Finally, consider the rest of the off-diagonal entries:

$\displaystyle(Z^{\prime\top}WZ^{\prime})_{ij}=Z^{\prime\top}_{i}WZ^{\prime}_{j}$
$\displaystyle=\sum_{S:i-1\in S,j-1\in S}w(S)$
$\displaystyle=\sum_{k=2}^{d}\binom{d-2}{k-2}\exp\left(-\left(1-\sqrt{k/d}\right)^{2}/\sigma^{2}\right)=a_{2},$

for $i,j=2,...,d+1$ such that $i\neq j$. Taken all these together,
$Z^{\prime\top}WZ^{\prime}$ has the closed form expression as shown in the
proposition.

We now proceed to find an expression for each element in $Z^{\prime\top}W$,
which has the first row entries:

$(Z^{\prime\top}W)_{1S}={\bm{1}}_{2^{d}}^{\top}W_{S}=w(S),$

for all $S\subseteq[d]$, whereas the rest of the rows have entries:

$(Z^{\prime\top}W)_{iS}=Z^{\prime\top}_{i}W_{S}=\mathbbm{1}\\{i-1\in
S\\}w(S),$

for $i=2,...,d+1$. Hence the closed form expression of $Z^{\prime\top}W$ in
the proposition follows. ∎

With closed form expressions for $Z^{\prime\top}WZ^{\prime}$ and
$Z^{\prime\top}W$, we numerically compute the LIME linear operator
$A_{\text{LIME}}=(Z^{\prime\top}WZ^{\prime})^{-1}_{[1:]}Z^{\prime\top}W$ and
compare it to the linear operators of the Banzhaf value ($A_{\text{Banzhaf}}$)
and Occlusion ($A_{\text{Occlusion}}$) with both large and small $\sigma^{2}$
values.

Figure 14: Frobenius norm of the difference between the LIME linear operator
$A_{\text{LIME}}$ and the Banzhaf linear operator $A_{\text{Banzhaf}}$, with
varying LIME exponential kernel hyperparameter $\sigma^{2}$ and number of
features. We include reference lines for the default LIME hyperparameters:
$\sigma^{2}=0.25$ for images, $\sigma^{2}=0.75\sqrt{d}$ for tabular data, and
$\sigma^{2}=25$ for text data.

As shown in Figure 14, as $\sigma^{2}$ increases, the Frobenius norm
difference between the LIME and Banzhaf linear operators approaches zero,
verifying our theoretical argument for the limiting case of
$\sigma^{2}\to\infty$. Furthermore, we observe that $A_{\text{LIME}}$ and
$A_{\text{Banzhaf}}$ have near zero difference for all default $\sigma^{2}$
values when there are $\\{10,15,20,25\\}$ features. As $\sigma^{2}$ approaches
zero, the Frobenius norm difference between the LIME and Occlusion linear
operators empirically approaches zero (Figure 15), confirming our intuition
for the limiting case where $\sigma^{2}\to 0$. Based on our computational
results, LIME with the default $\sigma^{2}$ values indeed has similar spectral
norm $\lVert\cdot\rVert_{2}$ and the operator norm
$\lVert\cdot\rVert_{1,\infty}$ compared to the Banzhaf value (Figure 16).

Figure 15: Frobenius norm of the difference between the LIME linear operator
$A_{\text{LIME}}$ and the Occlusion linear operator $A_{\text{Occlusion}}$,
with varying LIME parameter $\sigma^{2}$ and number of features. Figure 16:
Computed norms of summary technique linear operators. Note that all lines
except for RISE overlap on the left-hand-side showing the operator norm.

##  Appendix I Additional results for synthetic experiments

Figure 17: Occlusion, Banzhaf, RISE, and LIME attribution differences under
input perturbation with various perturbation norms in our synthetic data
experiment. Lines bounding the maximum attribution difference at each
perturbation norm are shown as empirical upper bounds. Theoretical upper
bounds are included for baseline and marginal feature removal, and omitted for
conditional feature removal because it is infinite for $\rho=1$.  Figure 18:
Slopes of empirical and theoretical upper bounds with respect to the input
perturbation norm, for Occlusion, Banzhaf, RISE, and LIME attribution
differences in the synthetic experiment. The parameter $\alpha$ varies the
logistic regression Lipschitz constant, and $\rho$ varies the correlation
between the two synthetic features ${\bm{x}}_{1},{\bm{x}}_{2}$. Figure 19:
Occlusion, Banzhaf, RISE, and LIME attribution differences between the
logistic regression classifiers $f,f^{\prime}$ as constructed in Section A.2
with baseline, marginal, and conditional feature removal, and varying
correlation $\rho$.

##  Appendix J Implementation details

Fully connected networks. For the wine quality dataset, we trained FCNs with
varying levels of weight decay. The architecture consists of $3$ hidden layers
of width $128$ with ReLU activations and was trained with Adam for $50$ epochs
with learning rate $0.001$. We also trained a conditional VAE (CVAE) following
the procedure in Frye et al. [27] to approximate conditional feature removal.
The CVAE encoder and decoder both consist of $2$ hidden layers of width $64$
with ReLU activations and latent dimension $8$. We trained the CVAE using a
prior regularization strength of $0.5$ for at most $500$ epochs with Adam and
learning rate $0.0001$, with early stopping if the validation mean squared
error did not improve for $100$ consecutive epochs. All models were trained
with NVIDIA GeForce RTX 2080 Ti GPUs with 11GB memory.

Convolutional networks. We trained CNNs for MNIST with varying levels of
weight decay. Following the LeNet-like architecture in Adebayo et al. [2], our
network has the following structure: Input $\rightarrow$ Conv($5\times 5$, 32)
$\rightarrow$ MaxPool($2\times 2$) $\rightarrow$ Conv($5\times 5$, 64)
$\rightarrow$ MaxPool($2\times 2$) $\rightarrow$ Flatten $\rightarrow$
Linear($1024,1024$) $\rightarrow$ Linear($1024,10$). We trained the model with
Adam for $200$ epochs with learning rate $0.001$. The same GPUs for FCN
training were used to train the CNNs.

ResNets. We trained ResNet-18 networks for CIFAR-10 from scratch with Adam for
$500$ epochs with learning rate $0.00005$. The same GPUs for FCN training were
used to train the ResNet-18 networks. We used a ResNet-50 network pre-trained
on ImageNet and fine-tuned it for Imagenette with Adam for $20$ epochs with
learning rate $0.0005$. The ResNet-50 networks were trained with NVIDIA Quadro
RTX6000 GPUs with 24GB memory.

##  Appendix K Additional results for practical implications

Figure 20: Test accuracy of FCN trained on the wine quality dataset and CNN
trained on MNIST with varying degree of weight decay. Figure 21: Test accuracy
of ResNet-18 trained on CIFAR-10 and ResNet-50 fine-tuned on Imagenette with
varying degree of weight decay. Figure 22: Attribution difference for networks
trained with increasing weight decay, under input perturbations with
perturbation norms $\\{0.1,0.5,1,2\\}$. The results include the wine quality
dataset with FCNs and baseline (replacing with training set minimums),
marginal, and conditional feature removal. Error bars show the mean and $95\%$
confidence intervals across explicand-perturbation pairs. Figure 23: Pearson
correlation of attributions for networks trained with increasing weight decay,
under input perturbations with perturbation norms $\\{0.1,0.5,1,2\\}$. The
results include the wine quality dataset with FCNs and baseline (replacing
with training set minimums), marginal, and conditional feature removal. Error
bars show the means and $95\%$ confidence intervals across explicand-
perturbation pairs. Figure 24: Spearman rank correlation of attributions for
networks trained with increasing weight decay, under input perturbations with
perturbation norms $\\{0.1,0.5,1,2\\}$. The results include the wine quality
dataset with FCNs and baseline (replacing with training set minimums),
marginal, and conditional feature removal. Error bars show the means and
$95\%$ confidence intervals across explicand-perturbation pairs. Figure 25:
Attribution difference for networks trained with increasing weight decay,
under input perturbations with perturbation norms $\\{0.1,0.5,1,2\\}$. The
results include MNIST with CNNs and baseline feature removal with either
training set means or zeros. Error bars show the means and $95\%$ confidence
intervals across explicand-perturbation pairs. Figure 26: Pearson correlation
of attributions for networks trained with increasing weight decay, under input
perturbations with perturbation norms $\\{0.1,0.5,1,2\\}$. The results include
MNIST with CNNs and baseline feature removal with either training set means or
zeros. Error bars show the means and $95\%$ confidence intervals across
explicand-perturbation pairs. Figure 27: Spearman rank correlation of
attributions for networks trained with increasing weight decay, under input
perturbations with perturbation norms $\\{0.1,0.5,1,2\\}$. The results include
MNIST with CNNs and baseline feature removal with either training set means or
zeros. Error bars show the means and $95\%$ confidence intervals across
explicand-perturbation pairs. Figure 28: Sanity checks for attributions using
cascading randomization for the CNN trained on MNIST. Features are removed by
replacing them with zeros. Attribution similarity is measured by Pearson
correlation and Spearman rank correlation. We show the means and $95\%$
confidence intervals across $10$ random seeds.  Figure 29: Sanity checks for
attributions using cascading randomization for the ResNet-18 trained on
CIFAR-10. Features are removed by replacing them with zeros. Attribution
similarity is measured by Pearson correlation and Spearman rank correlation.
We show the means and $95\%$ confidence intervals across $10$ random seeds.
Missing points correspond to undefined correlations due to constant
attributions under model randomization.


remarkRemark hypothesisHypothesis claimClaim

# Krylov Subspace Recycling With Randomized Sketching For Matrix Functions

Liam Burke School Of Mathematics, Trinity College Dublin, College Green,
Dublin 2, Ireland<EMAIL_ADDRESS>Stefan Güttel. S. G. acknowledges a
fellowship from the Alan Turing Institute under the EPSRC grant EP/W001381/1
and a Royal Society Industry Fellowship IF/R1/231032 Department of
Mathematics, The University of Manchester, M13 9PL Manchester, United Kingdom,
<EMAIL_ADDRESS>
###### Abstract

A Krylov subspace recycling method for the efficient evaluation of a sequence
of matrix functions acting on a set of vectors is developed. The method
improves over the recycling methods presented in [Burke et al.,
arXiv:2209.14163, 2022] in that it uses a closed-form expression for the
augmented FOM approximants and hence circumvents the use of numerical
quadrature. We further extend our method to use randomized sketching in order
to avoid the arithmetic cost of orthogonalizing a full Krylov basis, offering
an attractive solution to the fact that recycling algorithms built from
shifted augmented FOM cannot easily be restarted. The efficacy of the proposed
algorithms is demonstrated with numerical experiments.

###### keywords:

matrix function, Krylov method, subspace recycling, randomized sketching

65F60, 65F50, 65F10, 68W20

##  1 Introduction

This paper is concerned with the development of Krylov subspace recycling
algorithms for the efficient computation of a sequence of matrix function
applications of the form

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

with matrices $A^{(i)}\in\mathbb{C}^{N\times N}$ and vectors
$\mathbf{b}^{(i)}\in\mathbb{C}^{N}$. This is a common problem arising in a
variety of scientific computing applications. Common examples include
simulations of quantum chromodynamics (QCD) which often require the evaluation
of the matrix sign function $f(z)=\mbox{sign}(z)$ on a sequence of slowly
changing Dirac matrices [24, 14], and the time integration of stiff systems of
ordinary differential equations (ODEs) using exponential integrators, which
requires evaluation of the matrix exponential $f(z)=\mbox{exp}(z)$ on a
sequence of matrices and vectors which may change only slightly between time
steps [8, 25]. Perhaps the most well-known example occurs when $f(z)=z^{-1}$,
in which case (1) is equivalent to solving a sequence of linear systems of
equations; see e.g., [31] and references therein.

In the special case where the matrices $A^{(i)}=A$ remain fixed, the number of
different vectors $\mathbf{b}^{(i)}$ is relatively small and they are all
available simultaneously, block Krylov subspace methods [15, 16] may be used
for evaluating (1). The recycling algorithms we develop in this work are of
particular interest in the case when many vectors $\mathbf{b}^{(i)}$ are
available in sequence rather than simultaneously and/or the matrices $A^{(i)}$
are slowly changing. The idea is to update an augmentation space from one
problem to the next in order to improve convergence. We stress that recycling
is different from deflated restarting [27], the latter assuming that the
matrix $A=A^{(i)}$ and vector $\mathbf{b}=\mathbf{b}^{(i)}$ are fixed while
updating an augmentation space after each restart.

The _recycled FOM for functions of matrices_ presented in [7] (and therein
abbreviated as r(FOM)2) is the first Krylov subspace method to treat a
sequence of matrix function applications of the form (1) using subspace
recycling. The method is based on augmented FOM approximants to the solution
of the shifted linear systems appearing in an integral representation of
$f(A)\mathbf{b}$, and it uses numerical quadrature for the approximation of
the integral. For notational convenience we refer to this algorithm as _rFOM
(quad)_ throughout this paper. Unlike standard shifted FOM, the residuals of
augmented shifted FOM approximants are not necessarily collinear when the
shift varies, making it difficult to restart rFOM (quad). The rFOM (quad)
algorithm thus suffers from excessive storage and orthogonalization costs when
the matrices $A^{(i)}$ are non-Hermitian. This is in addition to the need of
deriving and evaluating an accurate quadrature formula for the integral
representation of $f(A)\mathbf{b}$.

The aim of this paper is to overcome the limitations of quadrature-based
recycling, developing a robust method that effectively makes use of the fact
that Krylov information may be reused across the different problems (1) to
speed up convergence.

The overall structure and key contributions are as follows.

  * •

In section 2 we derive new closed-form expressions for augmented subspace
approximations of matrix functions. The resulting method, referred to simply
as rFOM, is mathematically equivalent to the integral-based augmented FOM
approximants in [7], but it does not require quadrature for its implementation
and is observed to be numerically more robust for larger Arnoldi cycle lengths
(see Figure 1). We also prove a convergence result for rFOM when applied for
Stieltjes functions of Hermitian positive definite matrices.

  * •

In section 3 we show how to combine rFOM with randomized sketching (resulting
in srFOM), allowing for the use of nonorthogonal bases of the Krylov and
augmentation spaces, thereby reducing computational complexity. We show how to
utilize the sketched Rayleigh–Ritz approach from [28] to update the recycling
subspace.

  * •

We also state a closed-form expression for an augmented GMRES-type approximant
in section 3 and its sketched counterpart. In comparison to the sGMRES
approximant proposed in [20], our closed-form GMRES-type approximant does not
require numerical quadrature for its evaluation. (However, both approximants
are _not_ mathematically equivalent.)

  * •

Section 4 discusses various implementation details of our methods, in
particular, a SVD-based stabilization approach that demonstrably improves the
numerical robustness of srFOM, and an error estimator that can be used to
dynamically control the dimension of the Krylov space used for each problem.

  * •

Finally, section 5 contains numerical experiments to demonstrate the efficacy
of the new rFOM and srFOM in comparison to the state-of-the-art.

For a visual illustration of the algorithmic main contributions in this paper
we refer to Figure 1 and its caption. The 20 QCD test problems solved in this
plot are described in more detail in section 5.1.

Figure 1: Solving $20$ QCD problems $(A^{(i)})^{-1/2}\mathbf{b}^{(i)}$ with
slowly changing matrices and varying vectors. We show the relative error
obtained for each problem with standard FOM (without recycling), the rFOM
(quad) algorithm of [7], and the new algorithms rFOM and srFOM. The Arnoldi
cycle length is $m=100$ for all problems, and the recycling methods use a
recycling subspace of dimension of $k=30$. The rFOM (quad) algorithm is run
using $50$ and $100$ quadrature nodes, yielding the same convergence curves.
This indicates that no improvement can be achieved by further increasing the
number of quadrature nodes. Algorithm 1 uses a closed form of the recycled FOM
approximant and does not require quadrature. Algorithm 2 additionally utilizes
randomized sketching to reduce orthogonalization cost and further speed up the
computation at the expense of a slightly increased error.

##  2 Generalized subspace extraction for matrix functions

The rFOM (quad) algorithm of [7] is derived from augmented FOM approximants
$\mathbf{x}_{m}(\sigma)$ for the shifted linear systems

$(\sigma I-A)\mathbf{x}(\sigma)=\mathbf{b}$

appearing in the integral representation

(2) $f(A)\mathbf{b}=\int_{\Gamma}(\sigma
I-A)^{-1}\mathbf{b}\,\mathrm{d}\mu(\sigma)=\int_{\Gamma}\mathbf{x}(\sigma)\,\mathrm{d}\mu(\sigma)$

with some (possibly complex-valued) measure $\,\mathrm{d}\mu$ defined on a
contour $\Gamma\subset\mathbb{C}$. The augmented approximation
$\mathbf{f}_{m}\approx f(A)\mathbf{b}$ is obtained by replacing
$\mathbf{x}(\sigma)\to\mathbf{x}_{m}(\sigma)$ in the integral representation
and evaluating the integral using numerical quadrature. Here we show that it
is not necessary to use quadrature, and that an augmented FOM algorithm and
many other variants of subspace extraction can be implemented by using a
closed-form expression involving the function $f$.

Let the columns of two matrices $\textbf{V}_{m}$ and $\textbf{W}_{m}$ span a
search and constraint space in $\mathbb{C}^{N}$, respectively. Then, following
[7, eq. (5.5)] with some notational simplifications, we define the
_Bubnov–Galerkin approximant for the solution of $(\sigma
I-A)\mathbf{x}(\sigma)=\mathbf{b}$_ as

(3) $\mathbf{x}_{m}(\sigma):=\mathbf{V}_{m}\left[\mathbf{W}_{m}^{*}(\sigma
I-A)\mathbf{V}_{m}\right]^{-1}\mathbf{W}_{m}^{*}\mathbf{b}.$

Clearly, $\mathbf{x}_{m}(\sigma)\in\operatorname{span}(\mathbf{V}_{m})$ and it
is easy to verify that the residual

$\mathbf{r}_{m}(\sigma):=\mathbf{b}-(\sigma I-A)\mathbf{x}_{m}(\sigma)$

is orthogonal to $\operatorname{span}(\mathbf{W}_{m})$:

$\displaystyle\mathbf{W}_{m}^{*}\mathbf{r}_{m}(\sigma)$ $\displaystyle=$
$\displaystyle\mathbf{W}_{m}^{*}[\mathbf{b}-(\sigma
I-A)\mathbf{x}_{m}(\sigma)]$ $\displaystyle=$
$\displaystyle\mathbf{W}_{m}^{*}[\mathbf{b}-(\sigma
I-A)\mathbf{V}_{m}\left[\mathbf{W}_{m}^{*}(\sigma
I-A)\mathbf{V}_{m}\right]^{-1}\mathbf{W}_{m}^{*}\mathbf{b}]$ $\displaystyle=$
$\displaystyle\mathbf{W}_{m}^{*}\mathbf{b}-\left[\mathbf{W}_{m}^{*}(\sigma
I-A)\mathbf{V}_{m}\right]\left[\mathbf{W}_{m}^{*}(\sigma
I-A)\mathbf{V}_{m}\right]^{-1}\mathbf{W}_{m}^{*}\mathbf{b}$ $\displaystyle=$
$\displaystyle\boldsymbol{0}.$

We also refer to [22] for an in-depth discussion of such (oblique) projection
approaches for matrix functions.

We stress that the approximant (3) is very general. In particular, when
$\mathbf{V}_{m}$ spans the Krylov space
$\mathcal{K}_{m}(A,\mathbf{b})=\operatorname{span}\\{\mathbf{b},A\mathbf{b},\ldots,A^{m-1}\mathbf{b}\\}$
and $\mathbf{V}_{m}=\mathbf{W}_{m}$, then (3) is the approximant produced by
the _full orthorgonalization method_ (FOM) [35]. If $\mathbf{V}_{m}$ spans the
Krylov space $\mathcal{K}_{m}(A,\mathbf{b})$ and $\mathbf{W}_{m}=(\sigma
I-A)\mathbf{V}_{m}$, then (3) is the _generalized minimal residual_ (GMRES)
approximant [36]. This is easy to verify by writing
$\mathbf{x}_{m}(\sigma)=\mathbf{V}_{m}\mathbf{y}_{m}(\sigma)$ and then
checking that the condition
$\mathbf{W}_{m}^{*}\mathbf{r}_{m}(\sigma)=\boldsymbol{0}$ is equivalent to the
normal equations for $\|\mathbf{b}-(\sigma
I-A)\mathbf{V}_{m}\mathbf{y}_{m}(\sigma)\|_{2}\to\min_{\mathbf{y}_{m}(\sigma)}$.
However, it is usually impractical to have a constraint space $\mathbf{W}_{m}$
that depends on the shift $\sigma$, so we also define a _GMRES-type
approximant_ by setting $\mathbf{W}_{m}=A\mathbf{V}_{m}$ independent of
$\sigma$.

Assuming that $\mathbf{W}_{m}^{*}\mathbf{V}_{m}$ is nonsingular, (3) can be
rewritten as

$\displaystyle\mathbf{x}_{m}(\sigma)$ $\displaystyle=$
$\displaystyle\mathbf{V}_{m}\left[\sigma\mathbf{W}_{m}^{*}\mathbf{V}_{m}-\mathbf{W}_{m}^{*}A\mathbf{V}_{m}\right]^{-1}\mathbf{W}_{m}^{*}\mathbf{b}$
$\displaystyle=$ $\displaystyle\mathbf{V}_{m}\left[\left(\sigma
I-\mathbf{W}_{m}^{*}A\mathbf{V}_{m}(\mathbf{W}_{m}^{*}\mathbf{V}_{m})^{-1}\right)\mathbf{W}_{m}^{*}\mathbf{V}_{m}\right]^{-1}\mathbf{W}_{m}^{*}\mathbf{b}$
$\displaystyle=$
$\displaystyle\mathbf{V}_{m}(\mathbf{W}_{m}^{*}\mathbf{V}_{m})^{-1}\left[\sigma
I-\mathbf{W}_{m}^{*}A\mathbf{V}_{m}(\mathbf{W}_{m}^{*}\mathbf{V}_{m})^{-1}\right]^{-1}\mathbf{W}_{m}^{*}\mathbf{b}.$

Alternatively, and equivalently,

(4) $\displaystyle\mathbf{x}_{m}(\sigma)$ $\displaystyle=$
$\displaystyle\mathbf{V}_{m}\left[\sigma
I-(\mathbf{W}_{m}^{*}\mathbf{V}_{m})^{-1}\mathbf{W}_{m}^{*}A\mathbf{V}_{m}\right]^{-1}(\mathbf{W}_{m}^{*}\mathbf{V}_{m})^{-1}\mathbf{W}_{m}^{*}\mathbf{b},$

and we prefer to use (4) for now. Using the integral representation (2) and
replacing $\mathbf{x}(\sigma)$ by the approximant $\mathbf{x}_{m}(\sigma)$
defined in (3), we obtain

(5)
$\mathbf{f}_{m}:=\mathbf{V}_{m}f\left((\mathbf{W}_{m}^{*}\mathbf{V}_{m})^{-1}\mathbf{W}_{m}^{*}A\mathbf{V}_{m}\right)(\mathbf{W}_{m}^{*}\mathbf{V}_{m})^{-1}\mathbf{W}_{m}^{*}\mathbf{b}.$

We refer to $\mathbf{f}_{m}$ as the _Bubnov–Galerkin approximant of
$f(A)\mathbf{b}$_; see also [22]. When reading this formula, it might be
helpful to keep in mind that
$\mathbf{V}_{m}(\mathbf{W}_{m}^{*}\mathbf{V}_{m})^{-1}\mathbf{W}_{m}^{*}$ is
the matrix representation of the oblique projector onto the span of
$\mathbf{V}_{m}$ along $\mathbf{W}_{m}$. It is remarkable that (5) is a closed
formula and its evaluation does not require quadrature as opposed to the
integral representation used in [7].

Another attractive property of (5) is that it lends itself to randomized
sketching, and we will exploit this in section 3 below. Indeed, a similar
closed-form representation has been derived for the sketched FOM approximant
in [20], but (5) is more general as it also encompasses GMRES and GMRES-type
approximants of $f(A)\mathbf{b}$.

###  2.1 Closed-form augmented FOM approximation

We now assume that $\mathbf{V}_{m}=\mathbf{W}_{m}$ with linearly independent
(but not necessarily orthonormal) columns. We then have

$(\mathbf{W}_{m}^{*}\mathbf{V}_{m})^{-1}\mathbf{W}_{m}^{*}=(\mathbf{V}_{m}^{*}\mathbf{V}_{m})^{-1}\mathbf{V}_{m}^{*}=\mathbf{V}_{m}^{\dagger},$

and therefore (5) simplifies to

(6)
$\mathbf{f}_{m}=\mathbf{V}_{m}f\left(\mathbf{V}_{m}^{\dagger}A\mathbf{V}_{m}\right)\mathbf{V}_{m}^{\dagger}\mathbf{b}.$

If the columns of $\mathbf{V}_{m}$ span the Krylov space
$\mathcal{K}_{m}(A,\mathbf{b})$, this is the Rayleigh–Ritz representation of
the FOM approximant as it has been studied extensively in the literature. In
the context of subspace recycling, however, $\mathbf{V}_{m}$ may span an
arbitrary augmentation space which is not necessarily a Krylov space.

The practical evaluation of (6) requires the computation of
$\mathbf{V}_{m}^{\dagger}A\mathbf{V}_{m}$. This can either be done by using
sketching (as we will discuss in section 3), or alternatively by
orthogonalization of $\mathbf{V}_{m}$. Following the augmented FOM case
studied in [7], $\mathbf{V}_{m}=\mathbf{W}_{m}=[U,V_{m}]$ contains an
augmentation basis $U\in\mathbb{C}^{N\times k}$ (such as $k$ approximants to a
few relevant eigenvectors of $A$) and $V_{m}\in\mathbb{C}^{N\times m}$ is a
Krylov basis. After orthogonalization
$\mathbf{V}_{m}\leftarrow\texttt{qr\\_econ}(\mathbf{V}_{m})$ we can simply
evaluate (6) as
$\mathbf{f}_{m}=\mathbf{V}_{m}f\left(\mathbf{V}_{m}^{*}A\mathbf{V}_{m}\right)\mathbf{V}_{m}^{*}\mathbf{b}$.

Our closed-form recycled FOM is given in the Algorithm 1 below, and Figure 1
demonstrates its numerical behavior on a sequence of QCD test problems. We
will discuss this problem in greater detail in section 5, but the key message
is that closed-form FOM with recycling can yield significant convergence
acceleration for sequences of matrix function computations, and it can also be
more stable than the quadrature-based rFOM (quad) [7]. We believe that the
improved stability of Algorithm 1 results, at least partially, from the
explicit orthogonalization of the augmented basis $\mathbf{V}_{m}$, which is
in contrast to recycling variants that split
$\mathbf{x}_{m}(\sigma)=\mathbf{s}_{m}(\sigma)+\mathbf{t}_{m}(\sigma)$ and
keep $U$ and $V_{m}$ separate, and also different from variants which work
with $\mathbf{V}_{m}$ but do not orthogonalize its columns.

######  Remark 2.1.

In cases where the dimension $m$ is known a priori or kept fixed, it can be
beneficial to flip the order of $V_{m}$ and $U$, using
$\mathbf{V}_{m}=\mathbf{W}_{m}=[V_{m},U]$. This is because the augmentation
matrix $U$ typically has fewer columns than $V_{m}$ (i.e., $k\ll m$) and if
the matrix $V_{m}$ has been computed by the Arnoldi process we already have an
Arnoldi decomposition

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

with $V_{m}^{*}V_{m}=I_{m}$, $V_{m}^{*}\mathbf{v}_{m+1}=0$, and an upper-
Hessenberg matrix $H_{m}\in\mathbb{C}^{m\times m}$. Hence ensuring that
$\mathbf{V}_{m}=[V_{m},U]$ has fully orthonormal columns requires the
orthonormalization of only $k$ additional vectors. Further, for orthonormal
$\mathbf{V}_{m}$ we have

(8)
$\mathbf{V}_{m}^{*}A\mathbf{V}_{m}=\begin{bmatrix}V_{m}^{*}AV_{m}&V_{m}^{*}AU\\\
U^{*}AV_{m}&U^{*}AU\end{bmatrix}=\begin{bmatrix}H_{m}&V_{m}^{*}AU\\\
h_{m+1,m}(U^{*}\mathbf{v}_{m+1})\mathbf{e}_{m}^{T}&U^{*}AU\end{bmatrix}.$

Hence, when care is taken to keep the Krylov basis numerically orthonormal
(e.g., using the modified Gram–Schmidt process with reorthogonalization), it
is possible to reduce the number of matrix-vector products with $A$ in
computing $\mathbf{V}_{m}^{*}A\mathbf{V}_{m}$.

In the special case when the matrices in the problem sequence (1) remain
fixed, it is possible to update an Arnoldi-like decomposition for $AU$ and
further reduce the number of matrix-vector products in forming (8). This is
the case, for example, when the new augmentation subspace $U^{(i+1)}$ is
obtained as

$U^{(i+1)}:=\mathbf{V}_{m}^{(i)}X=[V_{m}^{(i)},U^{(i)}]X$

from a previous cycle with some $X\in\mathbb{C}^{(m+k)\times k}$. If the
product $Y^{(i)}:=AU^{(i)}$ is available from the previous cycle in which have
also computed an Arnoldi decomposition
$AV_{m}^{(i)}=V_{m}^{(i)}H_{m}^{(i)}+h_{m+1,m}^{(i)}\mathbf{v}_{m+1}^{(i)}\mathbf{e}_{m}^{T}$,
then

$AU^{(i+1)}=\left[V_{m}^{(i)}(H_{m}^{(i)}X)+h_{m+1,m}^{(i)}\mathbf{v}_{m+1}^{(i)}\mathbf{e}_{m}^{T}X\,,\,Y^{(i)}X\right],$

which requires no additional matrix-vector products with $A$.

Algorithm 1 Closed-form recycled FOM for $f({A}^{(i)})\mathbf{b}^{(i)}$ rFOM

1: Input: ${A}^{(i)}\in\mathbb{C}^{N\times N}$,
$\mathbf{b}^{(i)}\in\mathbb{C}^{N}$, integers $m$ and $k$

2: Output: Approximants $\widehat{\mathbf{f}}_{m}^{(i)}\approx
f({A}^{(i)})\mathbf{b}^{(i)}$

3: Initialize $U\leftarrow[\ ]$

4: for $i=1,2,\ldots$ do

5: Generate decomposition
$A^{(i)}V_{m}=V_{m}H_{m}+h_{m+1,m}\mathbf{v}_{m+1}\mathbf{e}_{m}^{T}$ for
$\mathcal{K}_{m}({A}^{(i)},\mathbf{b}^{(i)})$

6: Set $\mathbf{V}_{m}\leftarrow[U,V_{m}]$

7: Orthonormalize $\mathbf{V}_{m}$ and construct
$\mathbf{V}_{m}^{*}A^{(i)}\mathbf{V}_{m}$

8:
$\widehat{\mathbf{f}}_{m}^{(i)}\leftarrow\mathbf{V}_{m}f(\mathbf{V}_{m}^{*}A^{(i)}\mathbf{V}_{m})\mathbf{V}_{m}^{*}\mathbf{b}^{(i)}$

9: Update $U$ by extracting $k$-dimensional subspace from
$\operatorname{span}(\mathbf{V}_{m})$

10: end for

###  2.2 Convergence of augmented FOM

For the case that $A$ is a Hermitian positive definite matrix and $f$ obeys
the Stieltjes integral representation

$f(z)=\int_{0}^{+\infty}(z+t)^{-1}\,\mathrm{d}\mu(t)$

we can essentially follow the arguments in [12] to prove convergence of the
augmented FOM approximant $\mathbf{f}_{m}$ defined in (6) to $f(A)\mathbf{b}$.
We spell out the details below.

Let us assume that $\mathbf{V}_{m}=[V_{m},U]$ where $V_{m}$ spans
$\mathcal{K}_{m}(A,\mathbf{b})$ and $U$ spans the augmentation space. Let us
introduce $\mathbf{x}(t):=(A+tI)^{-1}\mathbf{b}$, $t\geq 0$, and the
associated augmented FOM approximant

$\mathbf{x}_{m}(t):=\mathbf{V}_{m}(\mathbf{V}_{m}^{\dagger}A\mathbf{V}_{m}+tI)^{-1}\mathbf{V}_{m}^{\dagger}\mathbf{b}$

as well as the non-augmented (standard) FOM approximant

$\underline{\mathbf{x}}_{m}(t):=V_{m}(V_{m}^{\dagger}AV_{m}+tI)^{-1}V_{m}^{\dagger}\mathbf{b}.\
\ $

(Note the bold $\mathbf{V}_{m}$ versus non-bold $V_{m}$ typesetting.)
Associated with $\mathbf{x}_{m}(t)$ are the residuals
$\mathbf{r}_{m}(t):=\mathbf{b}-(A+tI)\mathbf{x}_{m}(t)$ and errors
$\mathbf{e}_{m}(t):=\mathbf{x}(t)-\mathbf{x}_{m}(t)$. The vectors
$\underline{\mathbf{r}}_{m}(t)$ and $\underline{\mathbf{e}}_{m}(t)$ are
defined analogously.

Denoting by $\lambda_{\max}$ the largest eigenvalue of $A$ and defining the
norm $\|\mathbf{e}\|_{A}:=\sqrt{\mathbf{v}^{*}A\mathbf{v}}$, it is easily
verified that for $t\geq 0$ we have

$\|\mathbf{v}\|_{A}\leq\sqrt{\lambda_{\max}/(\lambda_{\max}+t)}\cdot\|\mathbf{v}\|_{A+tI}.$

Further, the orthogonal residual conditions on $\mathbf{x}_{m}(t)$ and
$\underline{\mathbf{x}}_{m}(t)$ are equivalent to the minimal error conditions

$\|\mathbf{e}_{m}(t)\|_{A+tI}\to\min_{\mathbf{x}_{m}(t)\in\operatorname{span}(\mathbf{V}_{m})}\quad\text{and}\quad\|\underline{\mathbf{e}}_{m}(t)\|_{A+tI}\to\min_{\underline{\mathbf{x}}_{m}(t)\in\operatorname{span}(V_{m})}.$

Hence,

$\|\mathbf{e}_{m}(t)\|_{A+tI}\leq\|\underline{\mathbf{e}}_{m}(t)\|_{A+tI}$

as the minimization criterion defining the left-hand side takes place over the
augmented subspace which contains the non-augmented Krylov space. Now, we are
ready to conclude

$\displaystyle\|f(A)\mathbf{b}-\mathbf{f}_{m}\|_{A}$ $\displaystyle=$
$\displaystyle\left\|\int_{0}^{+\infty}\mathbf{e}_{m}(t)\,\mathrm{d}\mu(t)\right\|_{A}\leq\int_{0}^{+\infty}\|\mathbf{e}_{m}(t)\|_{A}\,\mathrm{d}\mu(t)$
$\displaystyle\leq$
$\displaystyle\int_{0}^{+\infty}{\frac{\sqrt{\lambda_{\max}}}{\sqrt{\lambda_{\max}+t}}}\cdot\|\mathbf{e}_{m}(t)\|_{A+tI}\,\mathrm{d}\mu(t)$
$\displaystyle\leq$
$\displaystyle\int_{0}^{+\infty}{\frac{\sqrt{\lambda_{\max}}}{\sqrt{\lambda_{\max}+t}}}\cdot\|\underline{\mathbf{e}}_{m}(t)\|_{A+tI}\,\mathrm{d}\mu(t).$

The last expression is an upper bound involving the error
$\underline{\mathbf{e}}_{m}(t)$ of the non-augmented (standard) FOM
approximant $\underline{\mathbf{x}}_{m}(t)$. It is precisely the expression
which is further bounded in [12] and the Lemma 4.1 and Corollary 4.4 therein.
In summary, we obtain the following version of [12, Cor. 4.4].

######  Theorem 2.2.

Let $A\in\mathbb{C}^{N\times N}$ be Hermitian positive definite,
$\mathbf{b}\in\mathbb{C}^{N}$, and let
$f(z)=\int_{t_{0}}^{+\infty}(z+t)^{-1}\,\mathrm{d}\mu(t)$ be a Stieltjes
function ($t_{0}\geq 0)$. Let $\mathbf{f}_{m}$ be the augmented FOM
approximant as defined in (6) and with $\mathbf{V}_{m}=[V_{m},U]$, where the
columns of $V_{m}$ span $\mathcal{K}_{m}(A,\mathbf{b})$ and $U$ spans an
augmentation space. Further, let $\lambda_{\min}$ and $\lambda_{\max}$ denote
the smallest and largest eigenvalues of $A$, respectively, and define the
functions

$\kappa(t):=\frac{\lambda_{\max}+t}{\lambda_{\min}+t},\quad
c(t):=\frac{\sqrt{\kappa(t)}-1}{\sqrt{\kappa(t)}+1},\quad\text{and}\quad\alpha_{m}(t)=\frac{1}{\cosh(m\ln
c(t))}.$

Then the augmented FOM approximant satisfies

$\|f(A)\mathbf{b}-\mathbf{f}_{m}\|_{A}\leq C\alpha_{m}(t_{0})$

with a constant $C=\|\mathbf{b}\|_{2}\sqrt{\lambda_{\max}}\cdot
f(\sqrt{\lambda_{\min}\lambda_{\max}})$.

Note that this theorem guarantees exactly the same convergence as [12, Cor.
4.4] for the non-augmented (standard) FOM. Without adding further conditions
on the augmentation space, which will likely be difficult to verify in
practice, it appears to be unclear how to quantify any convergence
improvements due to the augmentation. For an informal explanation of why
convergence improvements can be expected when the augmentation space contains
(approximate) eigenvectors of $A$ so that spectral deflation takes place, we
refer to [11, Sec. 4].

###  2.3 Closed-form augmented GMRES-type approximation

As mentioned previously, a GMRES-type approximant is obtained from (5) if
$\mathbf{W}_{m}=A\mathbf{V}_{m}$. In the case where $\mathbf{V}_{m}=[V_{m}]$
is an orthonormal basis of the Krylov space $\mathcal{K}_{m}(A,\mathbf{b})$
satisfying an Arnoldi decomposition (7), (5) can be written as

$\displaystyle\mathbf{f}_{m}$ $\displaystyle=$ $\displaystyle
V_{m}f\left((V_{m}^{*}A^{*}V_{m})^{-1}V_{m}^{*}A^{*}AV_{m}\right)(V_{m}^{*}A^{*}V_{m})^{-1}(AV_{m})^{*}\mathbf{b}$
$\displaystyle=$ $\displaystyle
V_{m}f\left(H_{m}+(|h_{m+1,m}|^{2}H_{m}^{-*}\mathbf{e}_{m})\mathbf{e}_{m}^{T}\right)V_{m}^{*}\mathbf{b}.$

This expression is well known (see, e.g., [19, eq. (1.5)] and [22, eq. (2.5)])
and it has been used in [12] (albeit with a typo) to show convergence of the
(restarted) harmonic Arnoldi method for positive real matrices.

##  3 Randomized sketching

As the dimension $m$ of the Krylov basis $V_{m}$ grows large, the cost of its
orthogonalization may dominate. In order to overcome this drawback, we may
work with a basis that is only partially orthogonalized and then use
randomized sketching to deal with the non-orthogonality. The representation
(6) is particularly well suited for sketching; see, e.g., [26, 4, 5, 2, 28,
3]. Assume that we have a matrix $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})$,

(9)
$(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., [37, 38, 26]. Condition (9)
can equivalently be stated with the Euclidean inner product [37, 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 (9) with high probability.

Using the sketching operator $S$, we can replace all the implicit least
squares problems in (6) by solutions of their sketched counterpart, leading to

(10)
$\widehat{\mathbf{f}}_{m}:=\mathbf{V}_{m}f\left([S\mathbf{V}_{m}]^{\dagger}[SA\mathbf{V}_{m}]\right)[S\mathbf{V}_{m}]^{\dagger}S\mathbf{b}.$

This representation may be badly affected by ill-conditioning of the sketched
basis $S\mathbf{V}_{m}$ and it is usually beneficial to perform basis
whitening [33] by computing an economic QR decomposition
$S\mathbf{V}_{m}=Q_{m}R_{m}$ and using

(11)
$\widehat{\mathbf{f}}_{m}=\mathbf{V}_{m}\left(R_{m}^{-1}f\left(Q_{m}^{*}[SA\mathbf{V}_{m}]R_{m}^{-1}\right)Q_{m}^{*}S\mathbf{b}\right).$

It is interesting to note that (11) is formally the same as the sketched FOM
approximant labeled (sFOM”’) in [20]. However, in that work $\mathbf{V}_{m}$
was chosen as a basis of the Krylov space $\mathcal{K}_{m}(A,\mathbf{b})$,
while here we allow for a basis of an arbitrary subspace of $\mathbb{C}^{N}$.

######  Remark 3.1.

The general form of the approximant (3) lends itself to sketching also in the
case where $\mathbf{W}_{m}=A\mathbf{V}_{m}$, i.e., for the GMRES-type
approximant. A closed-form sketched GMRES-type approximant is

(12)
$\widetilde{\mathbf{f}}_{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{b}.$

Note that all the sketched matrices appearing within $f(\cdot)$ are small
(unrelated to the original problem dimension $N$) and this formula is indeed
relatively cheap to evaluate. It may hence be an attractive alternative to the
quadrature-based sGMRES approximant presented in [20]. However, care must be
taken with the numerical implementation as the condition number of the
inverted matrices can become large. We have only performed some preliminary
tests with the sketched GMRES-type approximant and empirically found it to
work well, but here we prefer not to discuss this approximant any further and
instead focus on sketched and recycled FOM.

###  3.1 Updating the augmentation space

We now turn our attention to the problem of computing a sequence of matrix
functions $f(A^{(i)})\mathbf{b}^{(i)}$ and discuss how to update the
augmentation space $U$ from one problem to the next. Say, for one problem
$f(A^{(i)})\mathbf{b}^{(i)}$ we have computed a Krylov basis $V_{m}$ of
$\mathcal{K}_{m}(A^{(i)},\mathbf{b}^{(i)})$ and we have an augmentation basis
$U\in\mathbb{C}^{N\times k}$ accumulated over previous problems. 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^{(i)}\mathbf{V}_{m}$ at our
disposal. To obtain an updated augmentation space $U$, we can now follow the
sRR approach presented in [28, Sec. 6.3] by solving the problem

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

A solution (of minimum Frobenius norm) is
$\widehat{M}=(S\mathbf{V}_{m})^{\dagger}(SA^{(i)}\mathbf{V}_{m})$ or,
exploiting the basis whitening factorization $S\mathbf{V}_{m}=Q_{m}R_{m}$ from
above,

$\widehat{M}=R_{m}^{-1}Q_{m}^{*}(SA^{(i)}\mathbf{V}_{m}).$

We can now compute a partial Schur decomposition

$\widehat{M}X=XT,\quad X\in\mathbb{C}^{(m+k)\times k},\quad
T\in\mathbb{C}^{k\times k},$

associated with appropriately selected eigenvalues of $\widehat{M}$. (If a
quasi-triangular real Schur form is used, it may be necessary to increase $k$
to $k+1$ to not tear apart any $2\times 2$ diagonal blocks of $T$.) Finally,
the augmentation space and its sketched version are updated as

$U\leftarrow\mathbf{V}_{m}X,\quad SU\leftarrow(S\mathbf{V}_{m})X.$

In addition, if the matrix $A^{(i)}=A^{(i+1)}$ remains unchanged for the next
problem, we can also update $SA^{(i+1)}U\leftarrow(SA^{(i)}\mathbf{V}_{m})X$
without computing any sketches or matrix-vector products with $A^{(i+1)}$. We
summarize the overall FOM procedure which combines sketching and recycling in
Algorithm 2.

Algorithm 2 Sketch-and-recycle FOM for $f({A}^{(i)})\mathbf{b}^{(i)}$ srFOM

1: Input: ${A}^{(i)}$, $\mathbf{b}^{(i)}$, integers $m,k$, sketch parameter
$s$ ($m<s\ll N$)

2: Output: Approximants $\widehat{\mathbf{f}}_{m}^{(i)}\approx
f({A}^{(i)})\mathbf{b}^{(i)}$

3: Draw a sketching matrix $S\in\mathbb{C}^{s\times N}$

4: Initialize $U\leftarrow[\ ]$

5: for $i=1,2,\dots$ do

6: Generate (non-orthogonal) basis $V_{m}$ of
$\mathcal{K}_{m}({A}^{(i)},\mathbf{b}^{(i)})$

7: Set $\mathbf{V}_{m}\leftarrow[V_{m},U]$

8: Obtain sketch $S\mathbf{V}_{m}$ (and $SA^{(i)}\mathbf{V}_{m}$ if $i=1$ or
$A^{(i)}\neq A^{(i-1)}$)

9: Compute thin QR factorization ${Q}_{m}{R}_{m}\leftarrow{S}\mathbf{V}_{m}$

10:
$\widehat{\mathbf{f}}_{m}^{(i)}\leftarrow\mathbf{V}_{m}\left({R}_{m}^{-1}f({Q}^{*}_{m}[{S}{A}^{(i)}\mathbf{V}_{m}]{R}_{m}^{-1}){Q}^{*}_{m}{S}\mathbf{b}^{(i)}\right)$

11: Compute $\widehat{M}\leftarrow
R_{m}^{-1}Q_{m}^{*}[SA^{(i)}\mathbf{V}_{m}]$

12: Compute partial Schur form $\widehat{M}X=XT$

13: Update $U\leftarrow\mathbf{V}_{m}X$, $SU\leftarrow S\mathbf{V}_{m}X$ (and
$SA^{(i)}U\leftarrow SA^{(i)}\mathbf{V}_{m}X$)

14: end for

##  4 Implementation

In this section we discuss a number of points relating to the practical
implementation of Algorithm 2.

###  4.1 Stabilization

If the matrix $S\mathbf{V}_{m}=S[V_{m},U]$ is poorly conditioned, then we can
employ the regularization proposed in [28]. Instead of computing the QR
factorization $S\mathbf{V}_{m}=Q_{m}R_{m}$, we compute an economic SVD
$S\mathbf{V}_{m}=L\Sigma J^{*}$. We define the truncated matrices
$L_{\ell}=L(\,:\,,1:\ell)$, $\Sigma_{\ell}=\Sigma(1:\ell,1:\ell)$,
$J_{\ell}=J(\,:\,,1:\ell)$, where $\ell$ is the largest integer such that the
singular values of $S\mathbf{V}_{m}$ satisfy

$\sigma_{1}\geq\sigma_{2}\geq\cdots\geq\sigma_{\ell}\geq\texttt{svdtol}$

for a given tolerance (like $\texttt{svdtol}=10^{-14}$). Replacing

$S\mathbf{V}_{m}\rightarrow L_{\ell},\quad SA\mathbf{V}_{m}\rightarrow
SA\mathbf{V}_{m}J_{\ell}\Sigma_{\ell}^{-1},\quad\mathbf{V}_{m}\rightarrow\mathbf{V}_{m}J_{\ell}\Sigma_{\ell}^{-1}\
\ \text{(only implicitly!)}$

in (10), we obtain the stabilized sketched FOM approximant

(14)
$\widehat{\mathbf{f}}_{m}^{\text{stab}}=\mathbf{V}_{m}\left(J_{\ell}\cdot\Sigma_{\ell}^{-1}f\big{(}L^{*}_{\ell}[SA\mathbf{V}_{m}]J_{\ell}\Sigma^{-1}_{\ell}\big{)}L^{*}_{\ell}(S\mathbf{b})\right).$

The recycling subspace can be updated by computing an ordered QZ decomposition
so that both of

$Q(L_{\ell}^{*}SA\mathbf{V}_{m}J_{\ell})Z,\quad Q\Sigma_{\ell}Z$

are upper-triangular matrices with the $k$ targeted eigenvalues appearing in
the upper-left block. (If a quasi-triangular real QZ form is used, it may be
necessary to increase $k$ to $k+1$ to not tear apart any $2\times 2$ diagonal
blocks.) The recycling subspace can then be updated as
$U\leftarrow\mathbf{V}_{m}J_{\ell}Z(\,:\,,1\,:\,k)$.

######  Remark 4.1.

The sketched GMRES-type approximant (12) could be stabilized in the same
manner, but it is more convenient to compute a truncated SVD of
$S\mathbf{W}_{m}=SA\mathbf{V}_{m}$ (instead of $S\mathbf{V}_{m}$):
$S\mathbf{W}_{m}\approx L_{\ell}\Sigma_{\ell}J_{\ell}^{*}$. Upon replacing

$S\mathbf{V}_{m}\rightarrow S\mathbf{V}_{m}J_{\ell}\Sigma_{\ell}^{-1},\quad
S\mathbf{W}_{m}\rightarrow
L_{\ell},\quad\mathbf{V}_{m}\rightarrow\mathbf{V}_{m}J_{\ell}\Sigma_{\ell}^{-1}\
\ \text{(only implicitly!)},$

the stabilized version of the approximant (12) takes the form

$\widetilde{\mathbf{f}}_{m}^{\text{stab}}=\mathbf{V}_{m}\left(J_{\ell}f\left([L_{\ell}^{*}S\mathbf{V}_{m}J_{\ell}]^{-1}\Sigma_{\ell}\right)[L_{\ell}^{*}S\mathbf{V}_{m}J_{\ell}]^{-1}L_{\ell}^{*}(S\mathbf{b})\right).$

###  4.2 Computational complexity

We now compare the cost of sketched-and-recycled FOM to sketched FOM for
solving a single problem in the sequence (1). Following the analysis in [20]
we assume that

  * •

an Arnoldi cycle length of $m$ and a recycling subspace dimension of $k$ is
used,

  * •

$V_{m}$ is computed using the $t$-truncated Arnoldi process where
$t=\mathcal{O}(1)$, and

  * •

the sketching parameter is chosen as $s=\mathcal{O}(m+k)$.

For non-Hermitian problems, the dominant arithmetic cost of the Arnoldi method
(aside the unavoidable matrix-vector products with $A$) are the
$\mathcal{O}(Nm^{2})$ arithmetic operations required to orthogonalize the
Krylov basis $V_{m}$. One of the main advantages of sFOM and srFOM is that
these methods can work with non-orthogonal Krylov bases, e.g., generated by
the $t$-truncated Arnoldi process. In that process, the $(j+1)$-st Krylov
basis vector is computed by projecting $\mathbf{w}_{j}:=A\mathbf{v}_{j+1}$
against the previous $t$ basis vectors

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

and then setting
$\mathbf{v}_{j+1}:=\widehat{\mathbf{w}}_{j}/\|\widehat{\mathbf{w}}_{j}\|$;
see, e.g., [34, Sec. 3.3]. Truncated Arnoldi requires only $\mathcal{O}(Nmt)$
arithmetic operations for the basis generation, i.e., the cost grows linearly
in the basis dimension $m$. In a recent work [21], the possibility of using
the sketched basis $SV_{m}$ (which is readily available) to select $t$ vectors
among _all_ previously computed vectors (not just the most recent ones) has
been explored. Other possibilities for constructing Krylov bases using a
limited number of inner products (or no inner products at all) exist as well,
e.g., based on recurrence relations for Chebyshev polynomials [23, Sec. 4] or
Newton polynomials [32, Sec. 4].

If the sketching operator $S$ is a subsampled random discrete cosine or
Fourier transform [39, 26], then the cost of sketching a matrix
$\mathbf{V}_{m}$ with $m+k$ columns is $\mathcal{O}(N(m+k)\mbox{log}(m+k))$.
If $\mathbf{V}_{m}=[U,V_{m}]$ and $V_{m}$ is constructed by the truncated
Arnoldi process, we can cheaply obtain the sketch $SAV_{m}$ using the relation

$SAV_{m}=(SV_{m})H_{m}+h_{m+1,m}(S\mathbf{v}_{m+1})\mathbf{e}_{m}^{T},$

so overall only $m$ vector sketches are needed. (This observation also applies
to sFOM [20] without recycling but has not been utilized in that work.)

Performing the thin QR factorization in line $9$ requires a total cost of
$\mathcal{O}(s(m+k)^{2})=\mathcal{O}((m+k)^{3})$, while for large enough $N$
forming the full approximant in line $10$ is dominated by the $O(N(m+k))$ cost
for the linear combination with the columns of $\mathbf{V}_{m}$. This
computation should hence be avoided whenever possible, working only with short
$(m+k)$-dimensional coefficient vectors
$\left({R}_{m}^{-1}f({Q}^{*}_{m}[{S}{A}\mathbf{V}_{m}]{R}_{m}^{-1}){Q}^{*}_{m}{S}\mathbf{b}\right)$
for tasks like error estimation (see the next subsection). The dominant cost
of updating the recycling subspace (line $13$) is in the
$\mathcal{O}(N(m+k)k)$ arithmetic operations required to construct $U$. The
full dominant cost of sketched-and-recycled FOM is thus
$\mathcal{O}(N(m+k)\log(m+k)+(m+k)^{3})$. When compared to the
$\mathcal{O}(Nm\log(m)+m^{3})$ arithmetic operations required in sketched FOM
[20], we see that the dominant cost of both methods is the same when we choose
an Arnoldi cycle length of $m^{\prime}=m-k$.

###  4.3 Error estimation

It is widely appreciated that the convergence analysis of Krylov methods for
nonsymmetric matrices is challenging. The difficulties will be even more
pronounced when recycling and sketching techniques are incorporated into the
algorithms. We do not currently have a practical error analysis of our
sketched-and-recycled FOM (or GMRES) algorithms. Even in the case of sketched
FOM or GMRES alone, existing error bounds are rarely useful in practice as
they tend to overestimate the actual error by a large margin and often require
information about the matrix (like the numerical range) that is not easily
available [20]. See also [30] for some interesting insights into sFOM.

Luckily, simple and practical a-posteriori error estimates are easy to derive,
such as the frequently used difference of two iterates, i.e.,

$\big{\|}f(A)\mathbf{b}-\widehat{\mathbf{f}}_{m}\big{\|}\approx\big{\|}\widehat{\mathbf{f}}_{m+d}-\widehat{\mathbf{f}}_{m}\big{\|}$

for a small integer $d\geq 1$. Following [20], it is possible to approximately
evaluate this stopping criterion without access to the full matrix
$\mathbf{V}_{m}$ and without forming $\widehat{\mathbf{f}}_{m+d}$ and
$\widehat{\mathbf{f}}_{m}$ explicitly. Using the assumption that $S$ is an
$\varepsilon$-subspace embedding for $\mathbf{V}_{m}$, satisfying (9), we have

(15)
${\frac{1}{\sqrt{1+\varepsilon}}\big{\|}S(\widehat{\mathbf{f}}_{m+d}-\widehat{\mathbf{f}}_{m})\big{\|}}\leq\big{\|}\widehat{\mathbf{f}}_{m+d}-\widehat{\mathbf{f}}_{m}\big{\|}\leq\frac{1}{\sqrt{1-\varepsilon}}\big{\|}S(\widehat{\mathbf{f}}_{m+d}-\widehat{\mathbf{f}}_{m})\big{\|}$

Writing $\widehat{\mathbf{f}}_{m}=\mathbf{V}_{m}\widehat{\mathbf{y}}_{m}$, we
obtain from (15) the bound

$\displaystyle\big{\|}\widehat{\mathbf{f}}_{m+d}-\widehat{\mathbf{f}}_{m}\big{\|}$
$\displaystyle\leq$
$\displaystyle\frac{1}{\sqrt{1-\varepsilon}}\big{\|}SV_{m+d}\widehat{\mathbf{y}}_{m+d}-SV_{m}\widehat{\mathbf{y}}_{m}\big{\|}$
$\displaystyle=$
$\displaystyle\frac{1}{\sqrt{1-\varepsilon}}\left\|SV_{m+d}\left(\widehat{\mathbf{y}}_{m+d}-\left[\begin{array}[]{c}\widehat{\mathbf{y}}_{m}\\\
\boldsymbol{0}_{d}\end{array}\right]\right)\right\|.$

For an estimator, the unknown embedding quantity $\varepsilon$ can be set to a
fixed generous constant (like $\varepsilon=0.99$), or it can be estimated by
keeping track of $\|S\mathbf{v}_{j}\|/\|\mathbf{v}_{j}\|$ as the Krylov basis
vectors $\mathbf{v}_{j}$ are generated ($j=1,\dots,m$).

##  5 Numerical experiments

In this section we present results of numerical experiments demonstrating the
effectiveness of the rFOM (Algorithm 1) and srFOM (Algorithm 2) introduced
here. We compare these methods to the standard Arnoldi approximation (FOM),
the quadrature-based rFOM (quad) [7, Algorithm 2], and sketched FOM (sFOM)
[20, Algorithm 1].

In all methods where recycling is used, the recycling subspace is updated with
$k$ approximate eigenvectors corresponding to eigenvalues closest to the
origin. In all cases where randomized sketching is employed, the sketching
operator is a subsampled randomized discrete cosine transform (see, e.g., [26,
Sec. 9.3]) and the Arnoldi truncation parameter is $t=2$. The sFOM algorithm
uses basis whitening, but we run it _without any stabilization_ to provoke
potential instabilities. Stabilization _is_ used in srFOM (stab) with a fixed
tolerance of $\texttt{svdtol}=10^{-14}$.

All experiments were performed in MATLAB111 Code available at
https://github.com/burkel8/srFOM R2023a on a Windows 11 HP laptop with an 11th
Gen Intel(R) Core processor with 2.80 GHz and 8 GB of RAM. All reported
runtimes are averages over 10 repetitions of an experiment.

###  5.1 Inverse square root

In this set of experiments we consider the function $f(z)=z^{-1/2}$. This
function arises in several applications, including Dirichlet-to-Neumann maps
(see, e.g., [10]) and quantum chromodynamics (QCD) (see, e.g., [24, 7]).
Numerical simulations of QCD and the computation of spectral projectors often
require the evaluation of (4) involving the matrix sign function, which can be
represented as $\operatorname{sign}(A)=f(A^{2})A\mathbf{b}$.

We first approximate a sequence of $20$ vectors of the form (1) where the
first matrix $A^{(1)}$ in the sequence is a complex non-Hermitian lattice QCD
matrix of size $N=3,072$ obtained from the SuiteSparse Matrix Collection [9]
(matrix ID $1597$), plus a shift of $6.0777\cdot I$ to move all eigenvalues
away from the negative real axis (the left-most real eigenvalue of the
unshifted matrix is $-6.0776$). The other matrices in the sequence take the
form

$A^{(i+1)}=A^{(i)}+10^{-8}M^{(i)},\quad i=1,2,\ldots,19,$

where each $M^{(i)}$ is a unit Gaussian random matrix with the same sparsity
pattern as $A^{(1)}$, and the vectors $\mathbf{b}^{(i)}$ are randomly
generated with unit Gaussian entries.

Let us first discuss Figure 1 from the introduction in some more detail. This
figure shows the relative error

$\texttt{relerr}=\frac{\|f(A^{(i)})\mathbf{b}^{(i)}-\mathbf{f}_{m}^{(i)}\|}{\|f(A^{(i)})\mathbf{b}^{(i)}\|}$

obtained after $m=100$ Arnoldi iterations of rFOM (Algorithm 1) and srFOM
(Algorithm 2). A recycling subspace dimension of $k=30$, sketching dimension
of $s=400$, and an Arnoldi truncation parameter $t=2$ are used. We also show
the relative error obtained from the standard FOM approximation and two runs
of rFOM (quad) [7] with $50$ and $100$ quadrature nodes, respectively. The
error curves demonstrate the effectiveness of rFOM and srFOM in reducing the
relative error as the sequence of problems progresses and the recycling
subspaces improves. Additionally, Figure 1 demonstrates that rFOM is more
numerically robust than rFOM (quad), particularly for larger values of $m$. We
believe this improved robustness is at least partially attributable to the
explicit orthogonalization of $\mathbf{V}_{m}$. Thus for the remainder of
these experiments, we no longer include rFOM (quad) in the comparisons.

We now repeat the same experiment for Figure 2, but also include the error
curves for sFOM and the stabilized srFOM. We see that the sFOM approximants
attain an error that is only slightly larger than that of FOM, and that
stabilization indeed helps to further reduce the error of srFOM.

Figure 2: Relative errors when approximating a sequence of $20$ vectors of the
form $(A^{(i)})^{-1/2}\mathbf{b}^{(i)}$ using FOM, sFOM, rFOM, srFOM, and
stabilized srFOM. The Arnoldi cycle length is fixed at $m=100$ for all 20
problems. The other parameters are $k=30$ (recycling subspace dimension),
$s=400$ (sketching dimension), and $t=2$ (Arnoldi truncation length).

In Figure 3 we now allow the number of Arnoldi iterations $m$ to vary for each
problem in the sequence. This is done by checking the relative error every
$d=10$ iterations and stopping when relerr is below $10^{-10}$. We then plot
the values of $m$ required for each problem. We first note the rather high
value of $m$ required by srFOM for the first problem and the irregular
behavior of sFOM for many of the problems. This is caused by the absence of
truncated SVD-based stabilization. (On the first problem where there is no
available augmentation space $U$, both sFOM and srFOM are equivalent.) The
rather stringent tolerance of $10^{-10}$ pushes the non-stabilized sketching
methods to the limit. If the tolerance is raised to $10^{-9}$ or stabilization
is employed, the irregular behavior goes away. This can be seen in Figure 3
for the srFOM (stab) method, which behaves much more regularly with $m$
steadily decreasing as the sequence of problems progresses. Indeed, the values
of $m$ required by srFOM (stab) are only marginally above those of rFOM (which
uses no sketching).

In Table 1 we record the total number of matrix-vector products with the
matrices $A^{(i)}$ (MAT-VEC’s), the number of inner products, the number of
vector sketches, and the total runtime. We see that sFOM and rFOM are able to
reduce the number of inner products and MAT-VEC’s when compared to FOM,
respectively, resulting in a moderate reduction of total runtime. However, a
much better overall performance is obtained with the _combination_ of
recycling and sketching (with stabilization). This is because srFOM (stab)
combines the best of both worlds: it reduces the number of MAT-VEC’s using
recycling and it reduces orthogonalization cost via sketching.

Figure 3:  Arnoldi cycle length $m$ needed to approximate a sequence of $20$
vectors $(A^{(i)})^{-1/2}\mathbf{b}^{(i)}$. Here the Krylov dimension $m$ is
chosen adaptively for each problem to reach a relative error below
$\texttt{reltol}=10^{-10}$. The error criterion is checked every $d=10$
iterations. The other parameters are $k=30$, $s=400$, and $t=2$. The rather
irregular behavior of sFOM (and srFOM on the first problem) is caused by the
lack of stabilization.  inverse square root, fixed reltol
---
|  FOM |  sFOM |  rFOM |  srFOM |  srFOM (stab)
MAT-VEC’s |  3,000 |  3,610 |  2,660 |  3,140 |  2,860
Inner Products |  229,500 |  10,810 |  942,285 |  7,690 |  6,850
Sketches |  0 |  3,630 |  0 |  3,160 |  2,880
Runtime (seconds) |  10.6 |  9.8 |  8.4 |  6.6 |  5.7
Table 1:  Computational resources required to approximate a sequence of $20$
vectors $(A^{(i)})^{-1/2}\mathbf{b}^{(i)}$ to reach a relative error below
$\texttt{reltol}=10^{-10}$. The error criterion is checked every $d=10$
iterations. The other parameters are $k=30$, $s=400$, and $t=2$.

###  5.2 Linear systems of equations

Many applications in scientific computing require the solution to a sequence
of slowly changing linear systems of the form

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

a special case of (1) with $f(z)=z^{-1}$. Such problems arise, for example,
with PDE-constraint optimization problems [6].

Although Krylov subspace recycling is a well established technique for solving
sequences of linear systems (see, e.g., [31, 6, 17]), srFOM appears to be the
first attempt to combine recycling with randomized sketching. Of course, an
efficient Krylov solver for linear systems should also incorporate some
preconditioning and potentially restarting techniques. We will leave this to
future work but still think it is worth demonstrating that srFOM is a good
starting point for developing efficient linear system solvers that combine
recycling and sketching.

In this experiment we solve a sequence of $30$ linear systems with randomly
generated right-hand sides $\mathbf{b}^{(i)}$ having unit Gaussian entries.
The matrix $A$ is the Neumann matrix in MATLAB’s matrix gallery of size
$N=10,609$, which we shift by $+0.001\cdot I$ to make it nonsingular.

In Figure 4 we plot the Arnoldi cycle length $m$ required for FOM, sFOM, rFOM,
srFOM and srFOM (stab) to solve each linear system in the sequence with a
relative error below $10^{-9}$. (To be consistent with the rest of the paper,
we use an error-based stopping criterion, not a residual-based one.) The
stopping criterion is monitored every $d=10$ iterations. All recycling methods
use a recycling subspace dimension of size $k=30$, and all sketching methods
use a sketching dimension of $s=900$ and an Arnoldi truncation parameter
$t=2$. It is clear from Figure 4 that rFOM and srFOM require significantly
fewer Arnoldi iterations to reach convergence than FOM and sFOM. Additionally,
Figure 4 again illustrates the importance of stabilization, with srFOM (stab)
requiring consistently fewer Arnoldi iterations than srFOM.

Figure 4:  Arnoldi cycle length $m$ needed to approximate a sequence of $30$
vectors of the form $(A^{(i)})^{-1}\mathbf{b}^{(i)},i=1,2,\ldots,30$. Here the
Krylov dimension $m$ is chosen adaptively for each problem to reach a relative
error below $\texttt{reltol}=10^{-9}$. The error criterion is checked every
$d=10$ iterations. The other parameters are $k=30$, $s=900$, and $t=2$.

In Table 2 we record the total number of matrix-vector products with the
matrix $A$ (MAT-VEC’s), the number of inner products, the number of vector
sketches, and the total runtime to solve the sequence of 30 problems. Note
that rFOM alone already leads to a significant reduction of MAT-VEC’s and
runtime over FOM, though the number of required inner products for the
orthogonalization is very large. This is mitigated when sketching is employed,
with srFOM (stab) again being the fasted method overall.

inverse function, fixed reltol
---
|  FOM |  sFOM |  rFOM |  srFOM |  srFOM (stab)
MAT-VEC’s |  13,420 |  13,970 |  5,510 |  9,580 |  7,140
Inner Products |  3,022,030 |  41,880 |  3,975,815 |  28,710 |  21,390
Sketches |  0 |  14,000 |  0 |  9,610 |  7,170
Runtime (seconds) |  84.2 |  42.7 |  29.8 |  25.8 |  22.0
Table 2:  Computational resources required to approximate a sequence of $30$
vectors $(A^{(i)})^{-1}\mathbf{b}^{(i)}$ to reach a relative error below
$\texttt{reltol}=10^{-9}$. The error criterion is checked every $d=10$
iterations. The other parameters are $k=30$, $s=900$, and $t=2$.

###  5.3 Exponential function

We now consider the function $f(z)=\exp(z)$ which arises, e.g., with the
exponential time integration of ordinary differential equations. While Krylov
recycling is most suited for functions with finite singularities, we
demonstrate that the methods introduced in this paper can be readily applied
to entire functions as well. We include this example also because the
exponential has been a particularly difficult function to work with in
quadrature-based restarting and recycling methods as there is no canonical
contour for the integral representation [13, 7].

For our test we use the FEM discretization of an advection–diffusion problem
described in [1, Section 6.3]. The matrix $A$ of size $N=2,157$ was generated
using the COMSOL Multiphysics software. We consider a sequence of 30 problems
$\exp(0.01\cdot A)\mathbf{b}^{(i)}$ where $\mathbf{b}^{(1)}$ is chosen at
random with unit Gaussian entries and every subsequent vector in the sequence
is chosen as $\mathbf{b}^{(i+1)}=\exp(0.01\cdot A)\mathbf{b}^{(i)}$. This
mimics an exponential time stepping method for the solution of
$\mathbf{u}^{\prime}(t)=A\mathbf{u}$, $\mathbf{u}(0)=\mathbf{b}^{(1)}$. The
dimension of the augmentation space is $k=50$, the sketching parameter is
$s=400$, and the Arnoldi truncation parameter is $t=2$.

In this test we use the error estimator described in section 4.3 to terminate
the Arnoldi iteration for each problem, checking the estimator every $d=10$
iterations for a relative error tolerance of $10^{-9}$. Figure 5 (left) shows
the number of Arnoldi iterations $m$ performed for each problem. Note that $m$
decreases for all methods as the problems progress, even for the non-recycled
FOM and sFOM. This is because the solution vectors $f(A)\mathbf{b}^{(i)}$ get
progressively easier to approximate for later time steps as certain
eigenvector components in the starting vector $\mathbf{b}^{(i)}$ get damped
out by the exponential $f(z)=\exp(0.01z)$. Nevertheless, the recycling methods
still achieve a considerable further reduction in $m$. On the right of Figure
5 we show the actual relative error of the computed approximants (we have used
MATLAB’s expm to compute the exact solutions, which is still feasible for this
problem size). We find that these errors are indeed below the target tolerance
of $10^{-9}$ as intended, even though the error estimator may lead to a larger
$m$ than necessary in some cases.

Figure 5: Arnoldi cycle length $m$ (left) and relative errors for
approximating a sequence of $30$ vectors $\exp(0.01\cdot A)\mathbf{b}^{(i)}$
with an adaptive control of $m$ using an error estimator with delay parameter
$d=10$ and a target tolerance of $10^{-9}$. The other parameters are $k=50$,
$s=400$, and $t=2$.

##  6 Conclusions and future work

We have presented a new Krylov subspace recycling method for matrix functions
(rFOM) which, in contrast to existing methods, does not require any numerical
quadrature and performs more numerically robust. With srFOM we have shown how
randomized sketching can easily be incorporated to further reduce the
orthogonalization costs. We have presented results of numerical experiments
which suggest that rFOM and srFOM can approximate a sequence of matrix
functions with improved accuracy over FOM or sFOM, requiring less
computational work and runtime.

One of the observations we have made but not discussed in the numerical
experiments with recycling is that in some cases it is not actually necessary
to form $A^{(i+1)}U$ nor $SA^{(i+1)}U$ explicitly without any delay in
convergence even when $A^{(i+1)}\neq A^{(i)}$. In such cases the number of
matrix-vector products per problem can be further reduced by $k$, the
dimension of the augmentation space. With such a modification, the sRR step
(13) for problem $i+1$ in the sequence becomes

$\widetilde{M}:=\arg\min_{M\in\mathbb{C}^{(m+k)\times(m+k)}}\|S[A^{(i)}U^{(i)},A^{(i+1)}V_{m}^{(i+1)}]-S[U^{(i)},V_{m}^{(i+1)}]M\|_{F},$

where $U^{(i)}$ is the augmentation space extracted for the $i$-th problem and
$V_{m}^{(i+1)}$ is a Krylov basis for
$\mathcal{K}_{m}(A^{(i+1)},\mathbf{b}^{(i+1)})$. This is an inexact
Rayleigh–Ritz procedure for the matrix $A^{(i+1)}$ if $A^{(i)}\neq A^{(i+1)}$,
but as the sketching with $S$ introduces inexactness anyway, one might get
away with it. The analysis of this is potential future work.

In future work we would like to explore how mixed-precision computations could
be used to further reduce the runtime of rFOM and srFOM. Mixed-precision can
be combined with recycling in multiple ways. For example, one might carry out
parts of the Arnoldi process for forming $V_{m}$ using different precisions as
in [2], one could only do the sketching in lower precision (similarly to [18]
which use a lower-precision sketch to obtain a preconditioner for LSQR), or
one could use mixed-precision iterative refinement in the context of linear
systems as in [29]. Also, we believe that a thorough comparison of the
proposed GMRES-type methods (with and without sketching) against more
established recycling methods like GCRO-DR [31] for sequences (shifted) linear
systems would be very interesting.

For the case of linear systems of equations, we believe that srFOM and srGMRES
can be easily combined with preconditioning. However, if the preconditioner is
already effective in reducing the number of required MAT-VEC’s and hence
orthogonalization cost significantly, then the advantage of randomized
sketching dimishes. Hence sketching might be particularly useful in situations
where an efficient preconditioner is not easily available; see also the
discussion in [28, Sec. 1.4]. We hope that such aspects will be explored in
future work.

## Acknowledgments

This work was jointly funded by an Irish Research Council Government Of
Ireland Postgraduate Scholarship and the Manchester Mathematical Sciences
(MiMS). We are grateful for discussions with Marcel Schweitzer and Kirk
Soodhalter.

## References

  * [1] M. Afanasjew, M. Eiermann, O. G. Ernst, and S. Güttel.  Implementation of a restarted Krylov subspace method for the evaluation of matrix functions.  Linear Algebra Appl., 429:2293–2314, 2008.
  * [2] O. Balabanov and L. Grigori.  Randomized block Gram-Schmidt process for solution of linear systems and eigenvalue problems.  Technical Report arXiv:2111.14641, 2021.
  * [3] O. Balabanov and L. Grigori.  Randomized Gram–Schmidt process with application to GMRES.  SIAM J. Sci. Comput., 44(3):A1450–A1474, 2022.
  * [4] 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.
  * [5] 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.
  * [6] M. Bolten, E. de Sturler, C. Hahn, and M. L. Parks.  Krylov subspace recycling for evolving structures.  Comput. Methods Appl. Mech. Engrg., 391:Paper No. 114222, 15, 2022\.
  * [7] L. Burke, A. Frommer, G. Ramirez-Hidalgo, and K. M. Soodhalter.  Krylov subspace recycling for matrix functions.  arXiv preprint arXiv:2209.14163, 2022.
  * [8] M. Caliari, F. Cassini, and F. Zivcovich.  BAMPHI: matrix-free and transpose-free action of linear combinations of $\varphi$-functions from exponential integrators.  J. Comput. Appl. Math., 423:114973, 2023.
  * [9] T. A. Davis and Y. Hu.  The University of Florida sparse matrix collection.  ACM TOMS, 38:1–25, 2011.
  * [10] V. Druskin, S. Güttel, and L. Knizhnerman.  Near-optimal perfectly matched layers for indefinite Helmholtz problems.  SIAM Review, 58(1):90–116, 2016.
  * [11] M. Eiermann, O. G. Ernst, and S. Güttel.  Deflated restarting for matrix functions.  SIAM J. Matrix Anal. Appl., 32(2):621–641, 2011.
  * [12] A. Frommer, S. Güttel, and M. Schweitzer.  Convergence of restarted Krylov subspace methods for Stieltjes functions of matrices.  SIAM J. Matrix Anal. Appl., 35(4):1602–1624, 2014.
  * [13] A. Frommer, S. Güttel, and M. Schweitzer.  Efficient and stable Arnoldi restarts for matrix functions based on quadrature.  SIAM J. Matrix Anal. Appl., 35:661–683, 2014.
  * [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, volume 15.  Springer Science & Business Media, 2000.
  * [15] A. Frommer, K. Lund, and D. B. Szyld.  Block Krylov subspace methods for functions of matrices.  Electron. Trans. Numer. Anal., 47:100–126, 2017.
  * [16] A. Frommer, K. Lund, and D. B. Szyld.  Block Krylov subspace methods for functions of matrices II: Modified block FOM.  SIAM J. Matrix Anal. Appl., 41(2):804–837, 2020.
  * [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] V. Georgiou, C. Boutsikas, P. Drineas, and H. Anzt.  A mixed precision randomized preconditioner for the LSQR solver on GPUs.  In International Conference on High Performance Computing, pages 164–181. Springer, 2023.
  * [19] S. Goossens and D. Roose.  Ritz and harmonic Ritz values and the convergence of FOM and GMRES.  Numer. Linear Algebra Appl., 6(4):281–293, 1999.
  * [20] 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.
  * [21] S. Güttel and I. Simunec.  A sketch-and-select Arnoldi process.  Technical Report arXiv:2306.03592, 2023.
  * [22] M. Hochbruck and M. E. Hochstenbach.  Subspace extraction for matrix functions.  Technical report, Case Western Reserve University, Department of Mathematics, Cleveland, 2005.
  * [23] W. D. Joubert and G. F. Carey.  Parallelizable restarted iterative methods for nonsymmetric linear systems. part I: Theory.  Int. J. Comput. Math., 44(1-4):243–267, 1992.
  * [24] F. Knechtli, M. Günther, and M. Peardon.  Lattice Quantum Chromodynamics: Practical Essentials.  Springer, 2017.
  * [25] Y. Y. Lu.  Computing a matrix function for exponential integrators.  J. Comput. Appl. Math., 161(1):203–216, 2003.
  * [26] P.-G. Martinsson and J. A. Tropp.  Randomized numerical linear algebra: Foundations and algorithms.  Acta Numer., 29:403–572, 2020.
  * [27] R. B. Morgan.  GMRES with deflated restarting.  SIAM J. Sci. Comput., 24(1):20–37, 2002.
  * [28] Y. Nakatsukasa and J. A. Tropp.  Fast & accurate randomized algorithms for linear systems and eigenvalue problems.  Technical Report arXiv:2111.00113, 2022.
  * [29] E. Oktay and E. Carson.  Mixed precision GMRES-based iterative refinement with recycling.  Technical Report arXiv:2201.09827, 2022.
  * [30] D. Palitta, M. Schweitzer, and V. Simoncini.  Sketched and truncated polynomial Krylov methods: Evaluation of matrix functions.  Technical Report arXiv:2306.06481, 2023.
  * [31] 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.
  * [32] B. Philippe and L. Reichel.  On the generation of Krylov subspace bases.  Appl. Numer. Math., 62(9):1171–1186, 2012.
  * [33] V. Rokhlin and M. Tygert.  A fast randomized algorithm for overdetermined linear least-squares regression.  Proc. Natl. Acad. Sci. USA, 105(36):13212–13217, 2008.
  * [34] Y. Saad.  Krylov subspace methods for solving large unsymmetric linear systems.  Math. Comput., 37(155):105–126, 1981.
  * [35] Y. Saad.  Iterative Methods for Sparse Linear Systems, 2nd edition.  SIAM, Philadelphia, 2000.
  * [36] 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.
  * [37] 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.
  * [38] D. P. Woodruff.  Sketching as a tool for numerical linear algebra.  Found. Trends Theor. Comput. Sci., 10(1–2):1–157, 2014.
  * [39] F. Woolfe, E. Liberty, V. Rokhlin, and M. Tygert.  A fast randomized algorithm for the approximation of matrices.  Appl. Comput. Harmon. Anal., 25:335–366, 2008.


# Physical and mathematical modeling in experimental papers: achieving
robustness of mathematical modeling studies

Vitaly V. Ganusov
Department of Microbiology, University of Tennessee, Knoxville, TN 37996, USA
Department of Mathematics, University of Tennessee, Knoxville, TN 37996, USA
National Institute for Mathematical and Biological synthesis, Knoxville, TN
37996, USA

###### Abstract

Development of several alternative mathematical models for the biological
system in question and discrimination between such models using experimental
data is the best way to robust conclusions. Models which challenge existing
theories are more valuable than models which support such theories.

Mathematical modeling has been especially strong in physical sciences where
theories, based on basic principles, were able to predict energies associated
with mass (famous equation $E=mc^{2}$) and existence of fundamental particles
(e.g., positron). There have been many applications of mathematical modeling
in biology but because it has been harder to come up with basic principles,
predictive successes of mathematical theory-based predictions in biology are
more modest than that in physics. With the development of novel experimental
methods to collect unique data spanning different levels of organization, from
gene sequences and intracellular processes in individual cells to populations
of organisms, there has been an impressive rise in the number of mathematical
modeling studies in biology. Bioinformatics, commonly defined as a science
dealing with analysis of sequencing data (but see (1)), and systems biology,
commonly defined as a science dealing with intracellular gene networks, are
two most well recognized advances driven by the availability of novel
datasets. Several recent reviews focus on the usefulness of mathematical
modeling for understanding complex biological systems (2, 3, 4, 5).
Specifically, Möbius & Laan (3) illustrated how mathematical modeling helps to
interpret experimental data using example of molecular motors. While Möbius &
Laan (3) do describe some pros and cons of mathematical modeling, there was
limited discussion on what makes a good mathematical modeling paper and when
mathematical modeling can lead to robust conclusions.

The authors did not cite and did not discuss an important philosophical paper
on the use of mathematical modeling in natural sciences, a paper that every
mathematical modeler should read, think about, and discuss with his/her peers
(6). In contrast with prevailing view in the mathematical modeling community
(e.g., (3)), Oreskes _et al._ (6) argued that for open natural systems (and
most if not all systems in biology are open) mathematical models cannot make
robust predictions. Moreover, based on philosophical definitions the authors
argued that verification of mathematical models is impossible and validation
of models, as commonly used in the literature, is nearly impossible (6). They
further argued that the use of mathematical models should be rather heuristic,
to understand limits of specific mechanisms when models are constrained by
experimental data. If most model predictions are not robust (7), how should
one proceed with best use of mathematical modeling?

One approach, which was not strongly emphasized by Möbius & Laan (3) is the
“method of multiple working hypothesis” by Chamberlin (8) augmented by the
“strong inference” of Platt (9). In essence, instead of developing a single
mathematical model for the phenomenon in question as it appears to be the norm
in the field, one must consider several alternative models and use
experimental data as a judge in discriminating between alternative hypotheses
(9, 10). In cases when data are not available for model discrimination,
analysis of alternative models can inform which experiments may need to be
performed to discriminate between alternative models. Furthermore, analysis of
multiple alternative models could also highlight limitations on predictions of
any individual model, thus, testing robustness of the model predictions which
for any given model a priori are expected to not be robust (6, 11).
Discrimination between alternative models using experimental data allows to
reject some of these models and rejection, rather than confirmation of models
is the strongest application of the scientific method (9, 12). Indeed, as
Oreskes _et al._ (6) highlighted “models which challenge existing theories are
more valuable than models which confirm them.” Confirmation of simplest
possible models is reasonable but showing that simple models do not work is
clearly more illuminating highlighting our gaps in understanding the system in
question (13).

As mathematical modeling is taking an integral part of biology it is time that
education and training of mathematical modelers is moving beyond development
and analysis of single models. Rather, development of multiple mathematical
models which are tested against experimental data and finding which
mathematical models are not able explain experimental data needs to become a
norm for the mathematical modeling publications in the 21st century.

## References

  * Hogeweg (2011) Hogeweg, P. 2011 The roots of bioinformatics in theoretical biology.  _PLoS Comput Biol_ , 7(3), e1002 021.  doi:10.1371/journal.pcbi.1002021.
  * Servedio _et al._ (2014) Servedio, M. R., Brandvain, Y., Dhole, S., Fitzpatrick, C. L., Goldberg, E. E., Stern, C. A., Cleve, J. V. & Yeh, D. J. 2014 Not just a theory–the utility of mathematical models in evolutionary biology.  _PLoS Biol_ , 12(12), e1002 017.  doi:10.1371/journal.pbio.1002017.
  * Möbius & Laan (2015) Möbius, W. & Laan, L. 2015 Physical and Mathematical Modeling in Experimental Papers.  _Cell_ , 163(7), 1577–1583.  doi:10.1016/j.cell.2015.12.006.
  * Torres & Santos (2015) Torres, N. V. & Santos, G. 2015 The (Mathematical) Modeling Process in Biosciences.  _Front Genet_ , 6, 354.  doi:10.3389/fgene.2015.00354.
  * Karr _et al._ (2015) Karr, J. R., Takahashi, K. & Funahashi, A. 2015 The principles of whole-cell modeling.  _Curr Opin Microbiol_ , 27, 18–24.  doi:10.1016/j.mib.2015.06.004.
  * Oreskes _et al._ (1994) Oreskes, N., Shrader-Frechette, K. & Belitz, K. 1994 Verification, validation, and confirmation of numerical models in the Earth sciences.  _Science_ , 263(5147), 641–646.  doi:10.1126/science.263.5147.641.
  * De Boer (2012) De Boer, R. J. 2012 Which of our modeling predictions are robust?  _PLoS Comput Biol_ , 8(7), e1002 593.  doi:10.1371/journal.pcbi.1002593.
  * Chamberlin (1890) Chamberlin, T. C. 1890 The Method of Multiple Working Hypotheses: With this method the dangers of parental affection for a favorite theory can be circumvented.  _Science_ , 15, 92–96.
  * Platt (1964) Platt, J. R. 1964 Strong Inference: Certain systematic methods of scientific thinking may produce much more rapid progress than others.  _Science_ , 146(3642), 347–353.  doi:10.1126/science.146.3642.347.
  * Hilborn & Mangel (1997) Hilborn, R. & Mangel, M. 1997 _The ecological detective: confronting models with data_.  Princeton University Press, Princeton NJ.
  * Kirk _et al._ (2015) Kirk, P. D. W., Babtie, A. C. & Stumpf, M. P. H. 2015 SYSTEMS BIOLOGY. Systems biology (un)certainties.  _Science_ , 350(6259), 386–388.  doi:10.1126/science.aac9505.
  * Popper (2002) Popper, K. 2002 _The Logic of Scientific Discovery_.  Routledge Classics. Taylor & Francis.
  * Noecker _et al._ (2015) Noecker, C., Schaefer, K., Zaccheo, K., Yang, Y., Day, J. & Ganusov, V. V. 2015 Simple mathematical models do not accurately predict early SIV dynamics.  _Viruses_ , 7, 1189–1217.