11institutetext<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
University of Florence, Italy

# DiffDefense: Defending against Adversarial Attacks via Diffusion Models
††thanks: This work was supported by the European Commission under European
Horizon 2020 Programme, grant number 951911—AI4Media

Hondamunige Prasanna Silva  Lorenzo Seidenari  Alberto Del Bimbo

###### Abstract

This paper presents a novel reconstruction method that leverages Diffusion
Models to protect machine learning classifiers against adversarial attacks,
all without requiring any modifications to the classifiers themselves. The
susceptibility of machine learning models to minor input perturbations renders
them vulnerable to adversarial attacks. While diffusion-based methods are
typically disregarded for adversarial defense due to their slow reverse
process, this paper demonstrates that our proposed method offers robustness
against adversarial threats while preserving clean accuracy, speed, and plug-
and-play compatibility.

Code at: https://github.com/HondamunigePrasannaSilva/DiffDefence.

###### Keywords:

Diffusion models Adversarial defense Adversarial Attack.

##  1 Introduction

Figure 1: Overview of our approach. Adversarial attacks happen in image space
by adding crafted noise to a pattern $x$, shifting classifier’s output to a
wrong class. DiffDefence starts by drawing a sample $z^{1}_{T}$, to diffuse
iteratively, for $T$ steps into a reconstruction
$x^{(i)}_{r}=z^{(i)}_{0}=G\left(z^{(i)}_{T}\right)$; we then optimize
$z^{(i)}_{T}$ so that the diffusion output for a given optimized pattern
$z_{T}^{(i+1)}$, lies closer to the original attacked sample. In the figure,
we drop the diffusion step subscript for readability purposes .

The susceptibility of machine learning models to adversarial attacks is a
major challenge in the field of artificial intelligence. While various
techniques have been proposed to enhance the robustness of classifiers against
such attacks, there is a pressing need for more effective and efficient
solutions. In recent years, generative models such as Generative Adversarial
Networks (GANs)[13], Diffusion Probabilistic Models[17] have emerged as a
promising approach to improve the resilience of machine learning models
against adversarial attacks.

Modern deep generative models have a common structural similarity: the
generation of novel patterns is usually performed by transforming some random
latent code $z$. Sampling $z\sim p(z)$, where $p(z)$ is a known distribution
(e.g. $\mathcal{N}(0,I)$ and then computing $G(z)$, where $G(\cdot)$ is a deep
neural network allows the generation of new data. Given a model $G(\cdot)$,
trained on clean data we can assume that attacked samples $x^{*}$ have a
different distribution, therefore finding a latent code $z^{*}$ able to
generate $x^{*}$ should be hard. Our approach builds on the idea that given
some attacked pattern $x^{*}=x+\epsilon$ where $x$ is an unknown clean pattern
and $\epsilon$ is a perturbation crafted to induce some classifier into a
mistake, we should be able to find some latent code $z^{*}$ for which
$G(z^{*})$ is closer to the unknown clean pattern $x$ than to the attacked one
$x^{*}$. In Fig. 2 it can be seen how an attack can add subtle patterns
(center) to a clean image(left) and how DiffDefense recover a correctly
classified example(right).

In this paper, we present a novel approach that leverages Diffusion Models to
enhance the resistance of machine learning classifiers to adversarial attacks.
Our proposed method involves reconstructing the input image using a reverse
process of a diffusion model (see Fig. 1 for details), which improves the
model’s ability to withstand adversarial attacks. We show that this approach
offers comparable speed and robustness to other generative model-based
solutions. Moreover, our proposed defense mechanism can be applied as a plug-
and-play tool to any classifier without compromising its accuracy, provided
that the diffusion model can generate high-quality images. Overall, our
approach holds promise as a viable alternative to other more complex to train
models, such as GANs, for defending against adversarial attacks on machine
learning models, owing to the benefits offered by Diffusion Models.

Our contribution is threefold:

  * •

We are the first to use recently successful Denoising Diffusion Probabilistic
models as a plug-in algorithm for reconstruction based adversarial defense.
Differently from [35, 26] our approach is based on reconstruction thus not
requiring backward and forward passes for each optimization step. Moreover,
DDPMs are more stable in training with respect to GANs which have also been
used as a reconstruction tool[30].

  * •

Thanks to a superior reconstructive and representational power, DDPMs require
less prototype embeddings and iterations to extract a clean pattern from the
attacked one, leading to higher inference efficiency with respect to [35, 26,
30].

  * •

Finally, our approach does not require to be trained on adversarial patterns
and can be used to detect attacked images.

Figure 2: Left original example from MNIST. Center result of the DeepFool
attack (wrongly classified). Right our DiffDefense reconstruction (correctly
classified).

##  2 Related Works

We now cover existing state of the art on modern generative models in
adversarial machine learning, forming the base of our reconstruction based
adversarial defense approach. We than discuss recent methods of adversarial
attacks and defense.

###  2.1 Generative Models

Generative models[19, 13, 17, 31] have emerged as a powerful class of machine
learning algorithms that can create new data samples with characteristics
similar to a given dataset. Their central idea is to learn the underlying
probability data distribution, which could then be used to generate new
patterns via sampling. Interestingly, these models have also proven to be
particularly effective in adversarial scenarios, showcasing their ability to
create samples that can attack classifiers[36, 34]. Aforementioned models have
also shown the capability to learn a semantically coherent latent embedding
space. This property has been exploited, in adversarial scenarios, to remove
attacks from patterns. By reconstructing[30, 20], purifying [26] the perturbed
sample and generating new samples to bolster adversarial training[16], these
models can significantly enhance the security of machine learning systems.
Recent research has highlighted the potential role of generative models in
adversarial defense, as their primary objective is to produce fake data that
closely resembles real data. In this paper we want to investigate the use of
Diffusion Models to bolster the robustness of models against adversarial
attacks. For a more thorough coverage of modern Deep Generative Models we
refer the reader to [4].

###  2.2 Adversarial Attacks and Defense

An adversarial attack is a process that aims at altering a classifier input
pattern in order to get the classifier to output a wrong prediction. So given
an input $x_{i}$ and a corresponding label $y_{i}$, and a classifier
$\mathcal{C}$ an attack method will aim at obtaining a $x^{*}$ such that
$\mathcal{C}(x^{*})\neq y_{i}$. Given this definition any method completely
replacing pattern $x$ with a different pattern would suffice in making a
classifier mistake the label. For this reason, a constraint on the
perturbation of the attacked pattern $x^{*}$ is also required. Therefore the
attacked pattern $x^{*}$ must be close to the original one
$||x-x^{*}||<\epsilon$.

Adversarial attacks can work in white-box and black-box scenarios. In the
white-box scenario the classifier is known to the attacker. The full knowledge
of the classification method implies that, for example for a neural network
all weights are known and the attacking method can leverage this knowledge.
White-box approaches may exploit the computation of the gradient of the model
loss with respect to its parameters for a specific input, such as the Fast
Gradient Signed Method (FGSM)[14]. PGD [22, 21] improves over [14] by a
refined attack generation obtained by iterative Projected Gradient Descent
(PGD). Instead, Deepfool[23] attempts to find the closest decision boundary to
then perturb the input in that direction. Combination of attacks have been
also proposed in [11], combining parameter free versions of PGD with
SquareAttack[2]. EOT+PGD combines the concepts of Expectation over
Transformation (EOT)[3] and Projected Gradient Descent (PGD) to improve its
effectiveness. Elasticnet [6] exploits a combination of L1 and L2
regularisation terms to provide an optimal trade-off.

In the black-box scenario, the attacking method does not have access to the
classifier, which is the most realistic setting. When an attacker has no
knowledge of a classifier’s architecture and weights, they can employ a query-
based approach to perturb the input without relying on the gradient by
applying a perturbation on the input until the classifier changes its output.
Depending on the model’s feedback for a given query, an attack can be
classified as a score-based[7] or hard-label[5] attack. In the score-based
setting,the attacker exploits the model’s output probabilities of each
decision. Several attacks have been created with this approach, such as Square
Attack[2], which selects localized square-shaped updates at random positions.
Pixel Attack [32] show attacks are possible even with a single pixel
perturbation.

In the hard label-based approach, the attacker exploits the model’s final
decision output. Recently, SIGN-OPT[9] improved a previous work [8] using
fewer queries (20K) to attack, being faster than the previous and obtaining a
similar error rate of white-box attacks, but remaining much slower than the
white-box. However, a query-based approach may not be as efficient as a white-
box attack. Instead, a black-box attack can be employed by using the
transferability[27, 28] of perturbed images to attack the target model. The
attacker can use a substitute model, where they have full knowledge, to
generate adversarial images using white-box attacks that can then be used to
attack the target model. Further coverage of adversarial attack and defense
technique can be found in [1].

The issue of adversarial attacks in machine learning has prompted the
development of various defense mechanisms, we can roughly identify three main
approaches: adversarial training, adversarial reconstruction, and adversarial
purification. Adversarial training methods improve the robustness of a model
against adversarial attacks by augmenting its training set with adversarial
examples. Introduced first by [14], adversarial training has become one of the
most successful defense against adversarial attacks [21, 15, 29], adversarial
training can also be enhanced using generative models for data augmentation
[16]. The main limitation of this approach is that is mainly protecting
classifiers from methods used in the adversarial training.

Adversarial reconstruction approaches leverage the projection of patterns onto
a learned latent manifold to regenerate the original input from its
adversarial counterpart. Generative models [30, 20] are a natural choice to
learn such latent representation and to obtain clean reconstructions out of
attacked patterns. Other approaches leveraged super-resolution networks[24] or
trained a reconstruction network to minimize the perceptual loss between the
reconstruction of the attacked pattern and the clean image.

Adversarial purification techniques perform a filtering of attacked patterns
removing adversarial perturbations while preserving their original features.
Recently proposed denoising Diffusion Models have been used as a tool for
purification [37, 33, 35]

Reconstructions based on GANs[30] are effective and generalize versus unseen
threats. However, the instability of GANs during training remains a challenge.
Moreover, many source noise embeddings and multiple reconstruction iterations
are required to obtain effective defense. Defense-VAE is faster and as
effective as Defense-GAN. However, to obtain effective reconstructions [20]
the method is fine-tuned on attacked images making the approach less general
and more prone to fail on unseen threats. Purification via Diffusion
Models[26, 33, 35] exploits multiple forward/backward passes to obtain a
reliable defense which requires significant time to purify an image.

With respect to [20, 30] our approach exploits powerful Diffusion Models as a
reconstruction tool. Interestingly, our approach is more efficient then
Defense-GAN, requiring less iterations and source embeddings. Different
from[20] we do not require to train on adversarial examples to work as
Defense-VAE[20]. Current defense mechanisms exploiting Diffusion models are
less efficient, requiring as much as 5s on a V100 card[26], while our approach
runs in 0.28s on a TitanXP card.

##  3 Methodology

We propose a diffusion reconstruction method as a defense against adversarial
attacks. The underlying idea is that adversarial attacks seek to deceive a
deep neural network (DNN) by introducing a disturbance to the image while
preserving its semantic meaning. Hence, the adversarial image ought to be
situated close to the original, unperturbed image.

Our approach is based on the idea that it is possible to induce a Generative
model $G(\cdot)$ to produce a given image $x^{*}$ by minimizing the distance
in image space of the output pattern, getting $\hat{z}$ as the result of such
minimization

$\hat{z}=\arg\min_{z}||G(z)-x^{*}||$ (1)

Obtaining the reconstructed image $x_{r}=G(\hat{z}$). Having $G(\cdot)$ being
learned on a clean dataset, the main assumption is that output generated
obtained by solving Eq. 1 are closer to clean examples than corrupted ones.

In our case $G(z)$ is the result of a reverse diffusion process, each step of
which is given by:

$\centering
z_{t-1}=\frac{1}{\sqrt{\alpha_{t}}}-\left(z_{t}-\frac{1-\alpha_{t}}{\sqrt{1-\bar{\alpha_{t}}}}\epsilon_{\theta}\left(z_{t},t\right)\right)+\sigma_{t}n\@add@centering$
(2)

where $\epsilon_{\theta}$ is the U-Net noise prediction model,
$\alpha_{t}=1-\beta_{t},\bar{\alpha_{t}}=\prod_{s=1}^{t}\alpha_{s},\\{\beta_{t}\in(0,1)\\}_{t=1}^{T},\sigma_{t}=\sqrt{\beta_{t}}$.

Our goal is to have the diffusion reverse process create a clean image, that
is as close as possible to the attacked input. To this end we must obtain a
suitable noise vector $z_{k}$. Therefore, we start from a random noise sample
$z_{T}^{1}$, and we iteratively generate an image using the reverse process of
a diffusion model. We then optimize z to solve Eq. 1 as shown in Algorithm 1.
In its general form, the proposed algorithm may also be run for multiple
source embeddings, although we found that it only increase slightly the
accuracy.

Algorithm 1 DiffDefense Reconstruction Algorithm. As a loss
$\mathcal{L}(x_{r}^{(i)},x^{*})$ we used Mean Square Error (MSE). $T^{*}$ are
the diffusion steps and $L$ are the gradient descent iterations, both treated
as hyperparameters. $\Delta=0.1$ is a decay rate.

1:Given adversarial image $x^{*}$

2:$z^{1}_{T}\sim\mathcal{N}(0,\,I)$

3:for $i=1,2,\ldots,L$ do

4: for $t=T^{*},T^{*}-1,\ldots,0\ steps$ do

5: $n\sim\mathcal{N}(0,\,I)$

6:
$z^{(i)}_{t-1}=\frac{1}{\sqrt{\alpha_{t}}}-\left(z^{(i)}_{t}-\frac{1-\alpha_{t}}{\sqrt{1-\bar{\alpha_{t}}}}\epsilon_{\theta}\left(z^{(i)}_{t},t\right)\right)+\sigma_{t}n$

7: end for

8: $\eta^{i}=\eta^{i-1}\Delta^{\frac{1}{\left\lceil L*0.8\right\rceil}}$

9: $x_{r}^{(i)}=z_{0}^{(i)}$

10:
$z_{T}^{(i+1)}=z_{T}^{(i)}-\eta^{i}\nabla_{z}\mathcal{L}(x_{r}^{(i)},x^{*})$

11:end for

###  3.1 Implementation details

The noise prediction U-Net $\epsilon_{\theta}$ architecture consists of a
contracting path, bottleneck layer and an expansive path. The contracting path
involves repeating a block with layer normalization, 3x3 convolutions, and
SiLU activation, followed by downsampling with a stride of 2. The number of
feature channels is doubled at each of the three downsampling step. The
bottleneck layer consists of the same block of the contracting path repeated
three times. The expansive path starts by concatenating the corresponding
feature map from the contracting path with an upsampled input using transpose
convolution. It is then followed by a block with layer normalization, 3x3
convolutions, and SiLU activation. At each upsampling step, the number of
feature channels is halved. This process is also repeated three times. Both
the contracting and expansive paths include a time-embedding layer, which
consists of two linear layers with a SiLU activation in between. This time-
embedding layer is added at each block of the contracting and expansive paths.

We employ two classifiers: the attacked classifier A and the surrogate
classifier B. Classifier A is composed of two 5x5 convolutions with 64 output
channels and stride of 2 and 1, respectively, using ReLU activations. It is
then followed by a dropout layer (p = 0.25), a linear layer with 128 output
features using ReLU activation, another dropout layer (p = 0.5), and finally a
linear layer with 10 output features. We use classifier B to generate
adversarial samples for black box attacks. B consists of a dropout layer (p =
0.2), followed by three convolutions with respective filter sizes of 8x8, 6x6,
and 5x5 and strides of 2, 2, and 1, using ReLU activations. Afterward, another
dropout layer (p = 0.5) is applied, and the final layer is a linear layer with
10 output features.

The Diffusion Model and the classifier are trained using the same clean
dataset. Training the classifier on reconstructed data is unnecessary if the
diffusion model generates high-fidelity images resembling the originals.

To implement the adversarial attack used to evaluate DiffDefense, we used
adversarial robustness toolkit[25] and torchattacks[18].

##  4 Experiments

This section presents the experiments that evaluate the proposed method using
both black-box and white-box attacks. First, we evaluate performance against
three classic attacks in both settings[14, 22, 23]. In these experiments we
seek optimal values for the number of iterations for gradient descent $L$, the
embedding set size $R$, and the diffusion step $T^{*}$. Then, keeping these
hyperparameter fixed we test DiffDefense against unseen attacks [2, 11, 38,
6]. Finally, we use this method to detect adversarial samples. The experiments
are conducted on two different datasets, MNIST[12] and KMNIST[10].

Table 1:  Performance of DiffDefense on white box & black box attacks on MNIST
& KMNIST datset. We report accuracy for each attack with and without defense.
For Black blox attacks, adversarial images has been crafted using a substitute
classifier.

Dataset |  Attack |  Type |  Without defense |  With defense
---|---|---|---|---
MNIST |  No attack |  - |  $99.14\%$ |  $99.06\%$
DeepFool |  White box |  $0.95\%$ |  $98.16\%$
Black box |  $97.17\%$ |  $98.86\%$
PGD |  White box |  $5.81\%$ |  $95.94\%$
Black box |  $51.28\%$ |  $97.18\%$
FGSM |  White box |  $23.72\%$ |  $89,95\%$
Black box |  $15.81\%$ |  $91.28\%$
KMNIST |  No attack |  - |  $95.18\%$ |  $94.38\%$
DeepFool |  White box |  $2.93\%$ |  $93.92\%$
Black box |  $92.16\%$ |  $93.92\%$
PGD |  White box |  $26.83\%$ |  $84.85\%$
Black box |  $58.43\%$ |  $91.49\%$
FGSM |  White box |  $37\%$ |  $79.5\%$
Black box |  $49.96\%$ |  $88.53\%$

###  4.1 Result of white-box & black-box attack

We investigate DiffDefense ability to withstand both white-box and black-box
attacks. To this end, we subject it to three potent white-box attacks: FGSM,
PGD, and Deep Fool. Furthermore, we evaluate the performance of DiffDefense
against these same attacks in the black-box setting, where we generate
adversarial samples using an auxiliary classifier to attack the target
classifier. In Tab. 1 we can see how all attacks are pretty effective in both
settings except for DeepFool used as a black box method. In general using
DiffDefense we can always recover a correct classification for almost all
attacked examples.

Table 2:  Robustness of DiffDefense against unseen threats on MNIST dataset.
Adversarial training using adversarial sample crafted by FGSM attack
$\epsilon=0.3$.

Attack | Type | W/O Defense | W/ Defense | Adv. Training
---|---|---|---|---
FGSM $\epsilon=0.3$ | White Box | $23.72\%$ | $89.95\%$ | $98.02\%$
PGD $\epsilon=0.3$ | White Box | $5.81\%$ | $95.94\%$ | $79.59\%$
Deep Fool | White Box | $0.95\%$ | $98.16\%$ | $5.81\%$
EOT+PGD $\epsilon=0.3$ | White Box | $24.57\%$ | $96.46\%$ | $89.22\%$
Square Attack | Black Box | $43.31\%$ | $97.31\%$ | $93.09\%$
AutoAttack | White Box | $1.26\%$ | $88.09\%$ | $45.86\%$
Elastic Net | White Box | $0.75\%$ | $95.75\%$ | $0.62\%$

###  4.2 Defense against unseen threats

One of the significant limitations of existing adversarial training defense
methods is their inability to effectively address previously unseen threats.
DiffDefense does not require to observe adversarial patterns to work,
nonetheless previous experiments were performed seeking the optimal values for
hyperparameters $L,R,T^{*}$. In order to assess the robustness of our proposed
approach to such unseen attacks we conducted evaluations using four different
attack techniques: Square Attack[2], Auto Attack[11] , EOT+PGD[38] attack and
Elastic Net [6], without tuning hyperparameters. Results in Tab. 2 indicate
that our method is robust against all four of these previously unseen attack
methods. Here we also test the behavior of adversarial training with samples
produced by FGSM with $\epsilon=0.3$. Interestingly, DiffDefense obtains high
accuracy even in cases in which adversarial training, is not helping at
all[23, 11, 6].

###  4.3 Ablation Studies

To evaluate the effectiveness and speed of our proposed approach, we conducted
an analysis of the three main hyperparameters, iteration number $L$, embedding
set size $R$ and diffusion step $T^{*}$. We found that the proposed method
does not need the same amount of steps of the Diffusion Model but it converges
with less steps, as shown in Fig.3. Moreover, in a comparison with Defense-
GAN[30], the results of our experiments revealed that our method achieved
convergence with fewer iteration steps and a smaller embedding set, while also
requiring less time to converge than the GAN-based method. This was evident in
the results presented in Tab.3, which show the superiority of our proposed
approach over Defense-GAN. The metrics used for comparison include robust
accuracy, which measures the accuracy after applying the defense, and time,
which indicates the duration to reconstruct a single image.

Table 3: Comparison with Defense-GAN[30]. Fewer iterations (L) and smaller
embedding set (R) in DiffDefense lead to faster convergence and reduced time.
All tests made on MNIST using white-box FGSM attack ($\epsilon=0.3$) and the
same classifier as [30]

Method |  L |  R |  Time |  Robust Acc
---|---|---|---|---
Defense Gan[30] |  25 |  10 |  0.086 |  $79.98\%$
100 |  1 |  0.273 |  $50.11\%$
100 |  10 |  0.338 |  $89.11\%$
200 |  10 |  0.675 |  $91.55\%$
Ours |  5 |  1 |  0.280 |  $87.78\%$
5 |  5 |  0.280 |  $\textbf{89.95}\%$

$1$$5$$10$$15$$0$$10$$20$$30$$40$$50$$60$$70$$80$$90$$100$Diffusion model
stepsAccuracyFGSM [14]PGD[22]DF[23]EOT+PGD[38]AA[11]SA[2]EE[6]

Figure 3: Accuracy analysis of the classifier after DiffDefense has been
applied on different diffusion steps. Using L = 4 and R = 5.

###  4.4 Attack detection

Interestingly the results of our study indicate that non-perturbed images are
reconstructed with greater ease in comparison to those subjected to
adversarial attacks. This is expected since the diffusion model and the
classifier are trained on the same data, facilitating the diffusion in the
reverse process phase using an unperturbed image to an adversarial image. This
ease of reconstruction is reflected in significantly smaller reconstruction
errors after an equal number of iterations. These findings suggest that the
reconstruction error may serve as a potential indicator of the presence of an
attack. In Fig. 4 we show ROC curves varying the diffusion step for [6, 23].
For all other methods[22, 38, 11, 2] we get AUC $\in$ [.99, 1].

(a) Elastic net

(b) Deep fool

Figure 4: Attack detection ROC curves for DiffDefense. In our experiments
FGSM, PGD, EOT+PGD, AutoAttack, Square Attack yielded a AUC $\in$ [.99, 1].

##  5 Conclusion

We proposed DiffDefense, a novel method that uses Diffusion models for
reconstruction, enhancing classifier robustness against attacks. Empirical
evaluation demonstrated its efficacy, speed, and potential as an alternative
to GAN-based methods and adversarial purification methods based on diffusion
models. We also showed that our approach is effective against previously
unseen attacks, highlighting its robustness to new attacks. Additionally, we
illustrated the usefulness of reconstruction as a tool for adversarial
detection. Our findings suggest that Diffusion based adversarial defense by
reconstruction is a promising path toward developing secure AI systems. We
believe that future work may further improve our method by adopting better
solvers for more accurate and faster reconstruction.

## References

  * [1] Akhtar, N., Mian, A.: Threat of adversarial attacks on deep learning in computer vision: A survey. Ieee Access 6, 14410–14430 (2018)
  * [2] Andriushchenko, M., Croce, F., Flammarion, N., Hein, M.: Square attack: a query-efficient black-box adversarial attack via random search. In: Proc. of ECCV. pp. 484–501. Springer (2020)
  * [3] Athalye, A., Engstrom, L., Ilyas, A., Kwok, K.: Synthesizing robust adversarial examples. In: International conference on machine learning. pp. 284–293. PMLR (2018)
  * [4] Bond-Taylor, S., Leach, A., Long, Y., Willcocks, C.G.: Deep generative modelling: A comparative review of vaes, gans, normalizing flows, energy-based and autoregressive models. IEEE transactions on pattern analysis and machine intelligence (2021)
  * [5] Brendel, W., Rauber, J., Bethge, M.: Decision-based adversarial attacks: Reliable attacks against black-box machine learning models. arXiv preprint arXiv:1712.04248 (2017)
  * [6] Chen, P.Y., Sharma, Y., Zhang, H., Yi, J., Hsieh, C.J.: Ead: elastic-net attacks to deep neural networks via adversarial examples. In: Proceedings of the AAAI conference on artificial intelligence. vol. 32 (2018)
  * [7] Chen, P.Y., Zhang, H., Sharma, Y., Yi, J., Hsieh, C.J.: Zoo: Zeroth order optimization based black-box attacks to deep neural networks without training substitute models. In: Proceedings of the 10th ACM workshop on artificial intelligence and security. pp. 15–26 (2017)
  * [8] Cheng, M., Le, T., Chen, P.Y., Yi, J., Zhang, H., Hsieh, C.J.: Query-efficient hard-label black-box attack: An optimization-based approach. arXiv preprint arXiv:1807.04457 (2018)
  * [9] Cheng, M., Singh, S., Chen, P., Chen, P.Y., Liu, S., Hsieh, C.J.: Sign-opt: A query-efficient hard-label adversarial attack. arXiv preprint arXiv:1909.10773 (2019)
  * [10] Clanuwat, T., Bober-Irizar, M., Kitamoto, A., Lamb, A., Yamamoto, K., Ha, D.: Deep learning for classical japanese literature (2018)
  * [11] Croce, F., Hein, M.: Reliable evaluation of adversarial robustness with an ensemble of diverse parameter-free attacks. In: International conference on machine learning. PMLR (2020)
  * [12] Deng, L.: The mnist database of handwritten digit images for machine learning research. IEEE Signal Processing Magazine 29(6), 141–142 (2012)
  * [13] Goodfellow, I., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., Bengio, Y.: Generative adversarial networks. Communications of the ACM 63(11), 139–144 (2020)
  * [14] Goodfellow, I.J., Shlens, J., Szegedy, C.: Explaining and harnessing adversarial examples. preprint, arXiv (2014)
  * [15] Gowal, S., Qin, C., Uesato, J., Mann, T., Kohli, P.: Uncovering the limits of adversarial training against norm-bounded adversarial examples. arxiv (2020), preprint
  * [16] Gowal, S., Rebuffi, S.A., Wiles, O., Stimberg, F., Calian, D.A., Mann, T.A.: Improving robustness using generated data. Advances in Neural Information Processing Systems 34 (2021)
  * [17] Ho, J., Jain, A., Abbeel, P.: Denoising diffusion probabilistic models. Advances in Neural Information Processing Systems 33, 6840–6851 (2020)
  * [18] Kim, H.: Torchattacks: A pytorch repository for adversarial attacks. preprint, arXiv (2020)
  * [19] Kingma, D.P., Welling, M.: An introduction to variational autoencoders. Foundations and Trends® in Machine Learning 12(4), 307–392 (2019)
  * [20] Li, X., Ji., S.: Defense-vae: A fast and accurate defense against adversarial attacks. In: Learning, M., in Databases:, K.D. (eds.) International Workshops of ECML PKDD 2019, Würzburg, Germany, September 16–20, 2019, Proceedings, Part II. Springer, International Publishing (2020)
  * [21] Madry, A., Makelov, A., Schmidt, L., Tsipras, D., Vladu, A.: Towards deep learning models resistant to adversarial attacks. In: International Conference on Learning Representations (2018)
  * [22] Madry, A., Makelov, A., Schmidt, L., Tsipras, D., Vladu, A.: Towards deep learning models resistant to adversarial attacks. arXiv preprint arXiv:1706.06083 (2017)
  * [23] Moosavi-Dezfooli, S.M., Fawzi, A., Frossard, P.: Deepfool: a simple and accurate method to fool deep neural networks. In: Proceedings of the IEEE conference on computer vision and pattern recognition (2016)
  * [24] Mustafa, A., Khan, S.H., Hayat, M., Shen, J., Shao, L.: Image super-resolution as a defense against adversarial attacks. IEEE Transactions on Image Processing 29, 1711–1724 (2019)
  * [25] Nicolae, M.I., Sinn, M., Tran, M.N., Buesser, B., Rawat, A., Wistuba, M., Zantedeschi, V., Baracaldo, N., Chen, B., Ludwig, H., M. Molloy, I., Edwards, B.: Adversarial robustness toolbox v1. 0.0. Tech. rep. (2018)
  * [26] Nie, W., Guo, B., Huang, Y., Xiao, C., Vahdat, A., Anandkumar, A.: Diffusion models for adversarial purification. arXiv preprint arXiv:2205.07460 (2022)
  * [27] Papernot, N., McDaniel, P., Goodfellow, I.: Transferability in machine learning: from phenomena to black-box attacks using adversarial samples. arXiv preprint arXiv:1605.07277 (2016)
  * [28] Papernot, N., McDaniel, P., Goodfellow, I., Jha, S., Celik, Z.B., Swami, A.: Practical black-box attacks against machine learning. In: Proceedings of the 2017 ACM on Asia conference on computer and communications security. pp. 506–519 (2017)
  * [29] Rebuffi, S.A., Gowal, S., Calian, D.A., Stimberg, F.: Wiles,o., and mann, t (2021), fixing data augmentation to improve adversarial robustness. arXiv preprint
  * [30] Samangouei, P., Kabkab, M., Chellappa, R.: Defense-gan: Protecting classifiers against adversarial attacks using generative models. preprint, arXiv (2018)
  * [31] Song, J., Meng, C., Ermon, S.: Denoising diffusion implicit models. preprint, arXiv (2020)
  * [32] Su, J., Vargas, D.V., Sakurai, K.: One pixel attack for fooling deep neural networks. IEEE Transactions on Evolutionary Computation 23(5), 828–841 (2019)
  * [33] Wang, J., Lyu, Z., Lin, D., Dai, B., Fu, H.: Guided diffusion model for adversarial purification. arXiv preprint arXiv:2205.14969 (2022)
  * [34] Wang, X., He, K., Hopcroft, J.E.: At-gan: A generative attack model for adversarial transferring on generative adversarial nets. preprint 3, arXiv (2019)
  * [35] Wu, Q., Ye, H., Gu., Y.: Guided diffusion model for adversarial purification from random noise. preprint, arXiv (2022)
  * [36] Xiao, C., Li, B., Zhu, J.Y., He, W., Liu, M., Song, D.: Generating adversarial examples with adversarial networks. arXiv preprint arXiv:1801.02610 (2018)
  * [37] Yoon, J., Hwang, S.J., Lee, J.: Adversarial purification with score-based generative models. In: International Conference on Machine Learning. PMLR (2021)
  * [38] Zimmermann, R.S.: Comment on. preprint, Adv-BNN: Improved Adversarial Defense through Robust Bayesian Neural Network”.” arXiv (2019)


# Augmented unprojected Krylov subspace methods from an alternative view of an
existing framework

Liam Burke111Corresponding Author Email<EMAIL_ADDRESS>School Of Mathematics,
Trinity College Dublin, College Green, Dublin 2  Kirk M. Soodhalter  School Of
Mathematics, Trinity College Dublin, College Green, Dublin 2

###### Abstract

Augmented Krylov subspace methods aid in accelerating the convergence of a
standard Krylov subspace method by including additional vectors in the search
space. These methods are commonly used in High-Performance Computing (HPC)
applications to considerably reduce the number of matrix vector products
required when building a basis for the Krylov subspace. In a recent survey on
subspace recycling iterative methods [Soodhalter et al, GAMM-Mitt. 2020], a
framework was presented which describes a wide class of such augmented
methods. The framework describes these methods broadly in a two step process.
Step one involves solving a projected problem via a standard Krylov subspace,
or other projection method, and step two then performs an additional
projection into the augmentation subspace. In this work we show that the
projected problem one must solve in step one has an equivalent unprojected
formulation. We then show how this observation allows the framework to be
adapted to describe unprojected augmented Krylov subspace methods. We
demonstrate this with two examples. We first show one can recover the R3GMRES
algorithm from this view of the framework, and then we use the framework to
derive the first unprojected augmented Full Orthogonalization Method (FOM).
This method allows one to recycle information from previous system solves, and
we thus denote it as unproj rFOM (unprojected recycled FOM).

The corresponding author would like to thank The Hamilton Scholarship for
generously funding this work as part of his PhD. studies

##  1 Introduction

It is well known that iterative Krylov subspace methods are the most practical
and efficient means of approximating the solution to a linear system of
equations who’s coefficient matrix is large, sparse and never explicitly
stored in memory. In High-Performance Computing (HPC) applications, it is of
interest to accelerate the convergence of these methods and reduce the number
of matrix vector operations required while still maintaining approximation
accuracy within the specified level of tolerance. One approach is to use
augmented Krylov subspace methods wherein the solution is approximated over
the sum of an iteratively built $j$-dimensional Krylov subspace
$\mathcal{V}_{j}$ and some other fixed $k$-dimensional subspace $\mathcal{U}$
known as the augmentation subspace. The subspace $\mathcal{U}$ is carefully
chosen to contain useful information to aid in accelerating the convergence of
the iterative solver, and is commonly taken to be the subspace spanned by
approximate eigenvectors corresponding to approximations to the $k$ smallest
eigenvalues of the coefficient matrix. The technique of Krylov subspace
recycling is a special case of an augmented Krylov subspace methods wherein
the subspace $\mathcal{U}$ comes from the solution of a previous system solve.
There is a large body of literature dedicated to the study of augmentation and
recycling (see eg. [1, 2, 3, 4]). In particular, we refer the reader to the
recent survey [5] and references therein which provide a broad overview of the
topic.

In this survey, a framework was presented which describes almost all augmented
Krylov subspace methods broadly in a two step process. Step one involves
applying a standard Krylov subspace, or other projection method, to a
projected problem. While step two involves another projection into the
augmentation subspace. The framework from this survey builds on earlier works
[1, 3, 4] and we refer the reader to these works for an understanding of how
the framework developed.

Despite the original view of this framework, not all augmented Krylov subspace
methods work by building a Krylov subspace with a projected matrix. Some
methods still use a Krylov subspace built from the matrix in its original
unprojected form. Although unprojected methods using approximate eigenvectors
will in general perform less favourably than projected methods, there are
still reasons to consider using an unprojected method. For example, particular
hardware constraints may render the projector two expensive to apply. In
addition they have advantages when dealing with the solution of ill-posed
inverse problems. An example of an unprojected method designed for this
setting is the R3GMRES algorithm [6]. This will be an important example used
throughout this note and for further details we refer the reader to [6].

In [7], it has been noted that although unprojected methods are not typically
described within the framework of [5], the R3GMRES algorithm can be cast in
this framework, and indeed viewing this method within this framework has
potential to lead to implementation benefits. Similar to the GCRO-DR algorithm
[2], the R3GMRES algorithm was formulated in terms of minimizing the $2$ norm
of the residual over the augmented Krylov subspace. It was shown in [7] that
this full minimization problem can be decoupled into solving a smaller
minimization problem over the Krylov subspace itself, followed by making an
additional correction from the space $\mathcal{U}$. This was done using the
normal equations of the full minimization problem and computing a block LU
factorization to isolate an equation for the correction from both Krylov
subspace $\mathcal{V}_{j}$, and from the augmentation subspace $\mathcal{U}$.
For more details we refer the reader to [7, Theorem 4.1].

The main goal of this paper is to show that the projected problem one must
solve in step one of the framework of [1, 5] has an equivalent unprojected
formulation. We prove this by showing that in exact arithmetic the
approximation obtained at step $j$ of these two problems is the same. We then
show that by viewing this projected problem in its unprojected form, the
framework can be extended to describe a wider class of unprojected augmented
Krylov subspace methods. We also note that the difference in describing a
projected method and an unprojected method within this framework is entirely
dependent on how we formulate the problem in step one. The correction obtained
in step two always remains the same.

We aim to show how one can easily derive an unprojected augmented Krylov
subspace method from this view of the framework using two examples. First, we
show that for choices of constraint subspaces corresponding to residual
minimizing methods, the R3GMRES algorithm arises naturally from this
framework, albeit not in its original form. Instead it reduces to its
equivalent formulation as presented in [7, Theorem 4.1]. Using the equivalent
unprojected form of the problem in step one we show that the algorithm arises
naturally from this view of the framework, and can thus be derived more easily
than the approach taken in the proof of [7, Theorem 4.1]. Secondly, we derive
the first unprojected augmented Full Orthogonalization Method (FOM). The
method we derive is also compatible with Krylov subspace recycling and thus we
will denote this method as unproj rFOM (unprojected recycled FOM). Since this
is a new method, we verify its effectiveness using a MATLAB implementation.
Our implementation constructs a $\mathcal{U}$ from a previous linear system
using ritz vectors extracted from the Hessenberg matrix of the previous system
solve. We show how this $\mathcal{U}$ accelerates the convergence of a new
system.

The layout of this paper is as follows. In the next section we present an
overview of standard Krylov subspace methods and recall both the Full
orthogonalization method (FOM) [8] and Generalized minimum residual method
(GMRES) [9]. In section 3 we give an overview of the framework for augmented
Krylov subspace methods as presented in [5]. In Section 4 we prove the
equivalence of the projected problem to an unprojected problem. We then show
how one can recover the R3GMRES algorithm [6] from this framework. We then
introduce a unprojected recycled FOM algorithm and demonstrate its
effectiveness using a MATLAB experiment. We then present our conclusions in
section 5.

#### Notation

We will use calligraphic letters such as $\mathcal{V}$ denote subspaces and
corresponding bold case capital letters, such as V to denote matrices with
columns forming a basis for that subspace. We may also use a subscript to
specify the dimension of the subspace (eg $\textbf{V}_{j}$ and
$\mathcal{V}_{j}$). When discussing projection methods in terms of arbitrary
search and constraint spaces we follow the same notation as in [5] and denote
the search space as some calligraphic letter and the constraint space as the
same letter with a tilde (eg. $\mathcal{V}_{j}$ and
$\tilde{\mathcal{V}}_{j}$). This is to ensure we are assuming the subspaces
are not necessarily the same. We will denote
$\Phi_{\mathcal{V},\tilde{\mathcal{V}}^{\perp}}$ as being the projector onto
$\mathcal{V}$ along $\tilde{\mathcal{V}}^{\perp}$ given explicitly as
$\Phi_{\mathcal{V},\tilde{\mathcal{V}}^{\perp}}=\textbf{V}_{j}(\tilde{\textbf{V}}^{*}_{j}\textbf{V}_{j})^{-1}\tilde{\textbf{V}}^{*}_{j}$.
We use $\textbf{e}_{j}$ to denote a vector of zeros with a $1$ in the $j^{th}$
position. The dimension of this vector should be clear from context.

##  2 An overview of Krylov subspace methods

Consider the problem of approximating the solution to a linear system

$\textbf{A}\textbf{x}=\textbf{b}$ (2.0.1)

where $\textbf{A}\in\mathbb{C}^{n\times n}$ and $\textbf{b}\in\mathbb{C}^{n}$,
and the matrix A is assumed to be large and sparse. In this section we
describe a well known residual projection framework for solving such systems.
More specifically, we then move to discuss Krylov subspace methods. For this
section we refer the reader to [8, 9] as well as chapter’s 5 and 6 of [10].

We wish to project the problem of solving the linear system (2.0.1) onto the
iteratively built $j$ dimensional subspace $\mathcal{V}_{j}$ while imposing
that the residual be orthogonal to another subspace of same dimension, denoted
$\tilde{\mathcal{V}}_{j}$. We exploit knowledge of an initial approximation
$\textbf{x}_{0}$ and thus we seek an approximation from the affine subspace
$\textbf{x}_{j}\in\textbf{x}_{0}+\mathcal{V}_{j}$. This approximation can thus
be written as $\textbf{x}_{j}=\textbf{x}_{0}+\textbf{t}_{j}$. Since the
correction $\textbf{t}_{j}\in\mathcal{V}_{j}$, we have
$\textbf{t}_{j}=\textbf{V}_{j}\textbf{y}_{j}$ with
$\textbf{y}_{j}\in\mathbb{C}^{j}$.

In a practical implementation of a projection method one builds the basis
$\textbf{V}_{j}$. Computing the full approximation at step $j$ thus requires
one to determine $\textbf{y}_{j}$, which can be done using the residual
orthogonality constraint

$\textbf{r}_{j}=\textbf{b}-\textbf{A}(\textbf{x}_{0}+\textbf{t}_{j})\perp\tilde{\mathcal{V}}_{j}$
(2.0.2)

This constraint yields a procedure for obtaining the correction from
$\mathcal{V}_{j}$ and corresponding updated solution given by the equations

$\textbf{y}_{j}=(\tilde{\textbf{V}}_{j}^{*}\textbf{A}\textbf{V}_{j})^{-1}\tilde{\textbf{V}}_{j}^{*}\textbf{r}_{0},\hskip
28.45274pt\textbf{x}_{j}=\textbf{x}_{0}+\textbf{V}_{j}\textbf{y}_{j}$ (2.0.3)

The corresponding residual can be updated via

$\textbf{r}_{j}=\textbf{b}-\textbf{A}\textbf{x}_{j}$ (2.0.4)

A practical implementation of a projection method requires one to
appropriately chooses subspace $\mathcal{V}_{j}$ and
$\tilde{\mathcal{V}}_{j}$. In this work we focus on Krylov subspace methods
wherein the subspace $\mathcal{V}_{j}$ is chosen to be the $j^{th}$ order
Krylov subspace generated from the matrix A and vector b defined by

$\mathcal{K}_{j}(\textbf{A},\textbf{b})=\\{\textbf{b},\textbf{A}\textbf{b},\textbf{A}^{2}\textbf{b},..,\textbf{A}^{j-1}\textbf{b}\\}$
(2.0.5)

We build a basis for this Krylov subspace using the Arnoldi algorithm. This
algorithm constructs the vectors forming the Krylov basis and orthogonalizes
each new vector against the previous using the Gram-Schmidt process. The
algorithm produces a matrix $\textbf{V}_{j}$ with columns forming an
orthonormal basis for $\mathcal{K}_{j}(\textbf{A},\textbf{b})$ and a an upper
Hessenberg matrix $\bar{\textbf{H}}_{j}$ which stores the orthogonalization
coefficients. The algorithm is outlined below.

Result: $\textbf{V}_{m}\in\mathbb{C}^{n\times m}$,
$\bar{\textbf{H}}_{m}\in\mathbb{C}^{(m+1)\times m}$

1 Given $\textbf{A}\in\mathbb{C}^{n\times n}$ and
$\textbf{b}\in\mathbb{C}^{n}$, cycle length $m$
$\textbf{v}_{1}=\textbf{b}/\|\textbf{b}\|$
for  _$j=1,..,m$_ do

2 $\textbf{w}_{j}:=\textbf{A}\textbf{v}_{j}$
for  _$i=1,..,j$_ do

3 $h_{ij}=\textbf{v}_{i}^{*}\textbf{w}_{j};$
$\textbf{w}_{j}:=\textbf{w}_{j}-h_{ij}\textbf{v}_{i}$

4 end for

5 $h_{j+1,j}=\|\textbf{w}_{j}\|_{2}$
$\textbf{v}_{j+1}=\textbf{w}_{j}/h_{j+1,j}$

6 end for

Algorithm 1 Arnoldi Algorithm

We will use the notation $\textbf{H}_{j}\in\mathbb{C}^{j\times j}$ to denote
the matrix $\bar{\textbf{H}}_{j}$ with last row deleted. At the end of the $j$
cycle Arnoldi process we have that the Arnoldi relation is satisfied

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

where $h_{j+1,j}$ is the entry at row $j+1$ and column $j$ of
$\bar{\textbf{H}}_{j}$.

We now discuss two well known Krylov subspace methods, both differing by their
choice of constraint space $\tilde{\mathcal{V}}_{j}$. The Full
Orthogonalization Method (FOM) [8] takes
$\tilde{\mathcal{V}}_{j}=\mathcal{V}_{j}=\mathcal{K}_{j}(\textbf{A},\textbf{b})$.
With this choice we see that $\textbf{y}_{j}$ can be written as

$\textbf{y}_{j}=(\textbf{V}_{j}^{*}\textbf{A}\textbf{V}_{j})^{-1}\textbf{V}^{*}_{j}\textbf{r}_{0}=(\textbf{V}_{j}^{*}(\textbf{V}_{j}\textbf{H}_{j}+h_{j+1,j}\textbf{v}_{j+1}\textbf{e}_{j}^{T}))^{-1}\textbf{V}^{*}_{j}\textbf{r}_{0}=\textbf{H}^{-1}_{j}(\|\textbf{r}_{0}\|\textbf{e}_{1})$
(2.0.7)

And thus computing $\textbf{y}_{j}$ amounts to solving the small $j\times j$
linear system

$\textbf{H}_{j}\textbf{y}_{j}=\|\textbf{r}_{0}\|\textbf{e}_{1}$ (2.0.8)

Similarly, the GMRES method [9] takes
$\tilde{\mathcal{V}}_{j}=\textbf{A}\mathcal{K}_{j}(\textbf{A},\textbf{b})$
leading to the solution of a small least squares problem for $\textbf{y}_{j}$

$\textbf{y}_{j}=\operatorname*{argmin}_{\textbf{z}\in\mathbb{C}^{j}}\left\|\|\textbf{r}_{0}\|\textbf{e}_{1}-\bar{\textbf{H}}_{j}\textbf{z}\right\|$
(2.0.9)

Once we have computed $\textbf{y}_{j}$, the solution can be updated via
$\textbf{x}_{j}=\textbf{x}_{0}+\textbf{V}_{j}\textbf{y}_{j}$. The Arnoldi
relation then allows for a cheap update of the residual

$\textbf{r}_{j}=\textbf{b}-\textbf{A}\textbf{x}_{j}=\textbf{r}_{0}-\textbf{V}_{j+1}\bar{\textbf{H}}_{j}\textbf{y}_{j}$
(2.0.10)

Since the basis grows at each iteration of the Arnoldi process, a practical
implementation of FOM or GMRES may need to be restarted in order to limit
memory consumption. One chooses a cycle length $j$, runs $j$ steps of the
Arnoldi algorithm and updates the solution and corresponding residual. The
basis $\textbf{V}_{j}$ is discarded and the process is restarted with the new
updated solution approximation. We outline the structure of restarted FOM and
GMRES here.

Input: $\textbf{A}\in\mathbb{C}^{n\times n}$ ,
$\textbf{x}_{0},\textbf{b}\in\mathbb{C}^{n}$ , tolerance $\epsilon$, cycle
length $j$
Result: $\textbf{x}_{j}$ such that
$\|\textbf{b}-\textbf{A}\textbf{x}_{j}\|\leq\epsilon$

1 Algorithm
$\textbf{r}\leftarrow\textbf{b}-\textbf{A}\textbf{x}_{0}$
while  _$\|\textbf{r}\| >\epsilon$_ do

2 Build a basis for $\mathcal{K}_{j}(\textbf{A},\textbf{r})$ using Arnoldi to
generate $\textbf{V}_{j+1}$ and $\bar{\textbf{H}}_{j}$.
if  _FOM_ then

3 Compute $\textbf{y}_{j}=\textbf{H}_{j}^{-1}(\|\textbf{r}\|\textbf{e}_{1})$

4 end if

5 if  _GMRES_ then

6 Compute
$\textbf{y}_{j}=\operatorname*{argmin}_{\textbf{z}\in\mathbb{C}^{j}}\left\|\|\textbf{r}\|\textbf{e}_{1}-\bar{\textbf{H}}_{j}\textbf{z}\right\|$

7 end if

8 $\textbf{x}\leftarrow\textbf{x}+\textbf{V}_{j}\textbf{y}_{j}$
$\textbf{r}\leftarrow\textbf{r}-\textbf{V}_{j+1}\bar{\textbf{H}}_{j}\textbf{y}_{j}$

9 end while

Algorithm 2 FOM and GMRES

For the GMRES algorithm the minimization can be solved using Givens rotations
which allow one to cheaply monitor the residual without computing it
explicitly. We do not focus on this detail in this work and refer the reader
to the references above.

##  3 A framework for augmented Krylov subspace methods

In this section we describe the framework of augmented Krylov subspace methods
in its most general form as presented in [5]. An augmented Krylov subspace
method aims to project the problem of solving the linear system (2.0.1) onto
the sum of an iteratively built Krylov subspace $\mathcal{V}_{j}$ and some
fixed $k$ dimensional subspace $\mathcal{U}$ carefully chosen to aid in
accelerating the convergence of the Krylov subspace method. Following the same
notation as in [5], the framework can be derived as follows.

Firstly, we note that the full approximation at step $j$ can be written as the
sum of a correction from the space $\mathcal{U}$,
$\textbf{s}_{j}\in\mathcal{U}$ and a correction from the Krylov space
$\mathcal{V}_{j}$, $\textbf{t}_{j}\in\mathcal{V}_{j}$. Exploiting knowledge of
an initial approximation $\textbf{x}_{0}$, the full solution can be written as

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

with $\textbf{z}\in\mathbb{C}^{k}$ and $\textbf{y}_{j}\in\mathbb{C}^{j}$. We
impose the residual at step $j$ be orthogonal to the sum of two subspaces, not
necessarily the same as the search space denoted as $\tilde{\mathcal{U}}$ and
$\tilde{\mathcal{V}}_{j}$. Again, we assume $\tilde{\mathcal{U}}$ is fixed and
$\tilde{\mathcal{V}}_{j}$ is built iteratively. The residual constraint is
written explicitly as

$\textbf{r}_{j}=\textbf{b}-\textbf{A}(\textbf{x}_{0}+\textbf{U}\textbf{z}+\textbf{V}_{j}\textbf{y}_{j})\perp\tilde{\mathcal{U}}+\tilde{\mathcal{V}}_{j}$
(3.0.2)

Application of this full residual constraint leads to a linear system for z
and $\textbf{y}_{j}$ given as

$\begin{bmatrix}\tilde{\textbf{U}}^{*}\textbf{A}\textbf{U}&\tilde{\textbf{U}}^{*}\textbf{A}\textbf{V}_{j}\\\
\tilde{\textbf{V}}^{*}_{j}\textbf{A}\textbf{U}&\tilde{\textbf{V}}^{*}_{j}\textbf{A}\textbf{V}_{j}\end{bmatrix}\begin{bmatrix}\textbf{z}\\\
\textbf{y}_{j}\end{bmatrix}=\begin{bmatrix}\tilde{\textbf{U}}^{*}\textbf{r}_{0}\\\
\tilde{\textbf{V}}^{*}_{j}\textbf{r}_{0}\end{bmatrix}$ (3.0.3)

A block LU factorization of this linear system allows one to decouple an
equation for $\textbf{y}_{j}$ independent of $\textbf{z}_{j}$ represented by
the second row of the following linear system

$\displaystyle\begin{bmatrix}\tilde{\textbf{U}}^{*}\textbf{A}\textbf{U}&\tilde{\textbf{U}}^{*}\textbf{A}\textbf{V}_{j}\\\
0&\tilde{\textbf{V}}^{*}_{j}\textbf{A}\textbf{V}_{j}-\tilde{\textbf{V}}^{*}_{j}\textbf{A}\textbf{U}(\tilde{\textbf{U}}^{*}\textbf{A}\textbf{U})^{-1}(\tilde{\textbf{U}}^{*}\textbf{A}\textbf{V}_{j})\end{bmatrix}\begin{bmatrix}\textbf{z}\\\
\textbf{y}_{j}\end{bmatrix}=\begin{bmatrix}\tilde{\textbf{U}}^{*}\textbf{r}_{0}\\\
\tilde{\textbf{V}}^{*}_{j}\textbf{r}_{0}-\tilde{\textbf{V}}^{*}_{j}\textbf{A}\textbf{U}(\tilde{\textbf{U}}^{*}\textbf{A}\textbf{U})^{-1}\tilde{\textbf{U}}^{*}\textbf{r}_{0}\end{bmatrix}$

This allows one to solve for $\textbf{y}_{j}$ easily. We can then directly
compute z in terms of $\textbf{y}_{j}$ using the first equation which gives

$\textbf{z}=(\tilde{\textbf{U}}^{*}\textbf{A}\textbf{U})^{-1}\tilde{\textbf{U}}^{*}\textbf{r}_{0}-(\tilde{\textbf{U}}^{*}\textbf{A}\textbf{U})^{-1}(\tilde{\textbf{U}}^{*}\textbf{A}\textbf{V}_{j})\textbf{y}_{j}$
(3.0.4)

We now examine the equation for $\textbf{y}_{j}$ in more detail. The second
equation can be arranged to give

$\displaystyle\tilde{\textbf{V}}^{*}_{j}(\textbf{I}-\textbf{A}\textbf{U}(\tilde{\textbf{U}}^{*}\textbf{A}\textbf{U})^{-1}\tilde{\textbf{U}}^{*})\textbf{A}\textbf{V}_{j}\textbf{y}_{j}=\tilde{\textbf{V}}^{*}_{j}(\textbf{I}-\textbf{A}\textbf{U}(\tilde{\textbf{U}}^{*}\textbf{A}\textbf{U})^{-1}\tilde{\textbf{U}}^{*})\textbf{r}_{0}$

We can define
$\Phi_{\textbf{A}\mathcal{U},\tilde{\mathcal{U}}^{\perp}}:=\textbf{A}\textbf{U}(\tilde{\textbf{U}}^{*}\textbf{A}\textbf{U})^{-1}\tilde{\textbf{U}}^{*}$
as being the projector onto $\textbf{A}\mathcal{U}$ along
$\tilde{\mathcal{U}}^{\perp}$. We see that finding $\textbf{y}_{j}$ from the
above equation is equivalent to solving the projected problem

Find $\textbf{t}_{j}\in\mathcal{V}_{j}$ as an approximate solution to the
problem
$(\textbf{I}-\Phi_{\textbf{A}\mathcal{U},\tilde{\mathcal{U}}^{\perp}})\textbf{A}\textbf{t}=(\textbf{I}-\Phi_{\textbf{A}\mathcal{U},\tilde{\mathcal{U}}^{\perp}})\textbf{r}_{0}$
(3.0.5) subject to $\textbf{r}_{j}\perp\tilde{\mathcal{V}}_{j}$

The framework is described in terms of two main steps. Step one involves
constructing and solving the projected problem (3.0.5) for $\textbf{y}_{j}$
and step 2 uses this $\textbf{y}_{j}$ to determine the correction z. These two
steps then combine to yield the full solution update, which can be written as

$\displaystyle\textbf{x}_{j}=\textbf{x}_{0}+\textbf{U}(\tilde{\textbf{U}}^{*}\textbf{A}\textbf{U})^{-1}\tilde{\textbf{U}}^{*}\textbf{r}_{0}+\textbf{V}_{j}\textbf{y}_{j}-\textbf{U}(\tilde{\textbf{U}}^{*}\textbf{A}\textbf{U})^{-1}(\tilde{\textbf{U}}^{*}\textbf{A}\textbf{V}_{j})\textbf{y}_{j}$

The corresponding full residual update is then

$\displaystyle\textbf{r}_{j}=\textbf{r}_{0}-\textbf{A}\textbf{U}(\tilde{\textbf{U}}^{*}\textbf{A}\textbf{U})^{-1}\tilde{\textbf{U}}^{*}\textbf{r}_{0}-\textbf{A}\textbf{V}_{j}\textbf{y}_{j}+\textbf{A}\textbf{U}(\tilde{\textbf{U}}^{*}\textbf{A}\textbf{U})^{-1}\tilde{\textbf{U}}^{*}\textbf{A}\textbf{V}_{j}\textbf{y}_{j}$

We now outlined more details of these steps. We begin with an initial
approximation $\textbf{x}_{0}$ and corresponding residual $\textbf{r}_{0}$. We
now must construct and solve the projected problem (3.0.5). The right hand
side of this projected problem can be obtained by projecting the initial
residual

$\hat{\textbf{r}}=\textbf{r}_{0}-\textbf{A}\textbf{U}(\tilde{\textbf{U}}^{*}\textbf{A}\textbf{U})^{-1}\tilde{\textbf{U}}^{*}\textbf{r}_{0}=(\textbf{I}-\Phi_{\textbf{A}\mathcal{U},\tilde{\mathcal{U}}^{\perp}})\textbf{r}_{0}$
(3.0.6)

This residual projection corresponds to the first two terms in the full
residual update. It can also be seen that projecting the residual in this way
performs a corresponding change in the initial solution $\textbf{x}_{0}$ given
by the first two terms in the full solution update.

In order to solve the projected problem, a natural choice of search space
$\mathcal{V}_{j}$ is the Krylov subspace corresponding to the projected
problem

$\mathcal{V}_{j}=\mathcal{K}_{j}((\textbf{I}-\Phi_{\textbf{A}\mathcal{U},\tilde{\mathcal{U}}^{\perp}})\textbf{A},\hat{\textbf{r}})$
(3.0.7)

In order to build this Krylov subspace, one never has to explicitly construct
the operator
$(\textbf{I}-\Phi_{\textbf{A}\mathcal{U},\tilde{\mathcal{U}}^{\perp}})\textbf{A}$.
Instead we can apply the original matrix A to each Arnoldi vector and the
action of
$(\textbf{I}-\Phi_{\textbf{A}\mathcal{U},\tilde{\mathcal{U}}^{\perp}})$ can be
done implicitly throughout the Arnoldi process. For an orthogonal projection,
this is done using an additional orthogonalization of the Arnoldi vectors. The
coefficients of this orthogonalization are stored in the matrix
$\textbf{B}_{j}=(\tilde{\textbf{U}}^{*}\textbf{A}\textbf{U})^{-1}\tilde{\textbf{U}}^{*}\textbf{A}\textbf{V}_{j}$,
which is also built throughout the process. For an oblique projection it can
be done in a similar, numerically stable way using appropriate choices for U
and C through a procedure which closely resembles Gram Schmidt.

The full solution update is then

$\textbf{x}_{j}=\hat{\textbf{x}}+\textbf{V}_{j}\textbf{y}_{j}-\textbf{U}\textbf{B}_{j}\textbf{y}_{j}$
(3.0.8)

and the full residual update is then

$\displaystyle\textbf{r}_{j}=\hat{\textbf{r}}-\textbf{A}\textbf{V}_{j}\textbf{y}_{j}+\textbf{A}\textbf{U}\textbf{B}_{j}\textbf{y}_{j}$
(3.0.9)

This residual updated can be performed more cheaply using the Arnoldi relation
we obtain from building the basis for
$\mathcal{K}_{j}((\textbf{I}-\Phi_{\textbf{A}\mathcal{U},\tilde{\mathcal{U}}^{\perp}})\textbf{A},\hat{\textbf{r}})$.
This Arnoldi relation is given by

$\displaystyle(\textbf{I}-\Phi_{\textbf{A}\mathcal{U},\tilde{\mathcal{U}}^{\perp}})\textbf{A}\textbf{V}_{j}=\textbf{V}_{j+1}\bar{\textbf{H}}_{j}$
$\displaystyle\implies\textbf{A}\textbf{V}_{j}=\Phi_{\textbf{A}\mathcal{U},\tilde{\mathcal{U}}^{\perp}}\textbf{A}\textbf{V}_{j}+\textbf{V}_{j+1}\bar{\textbf{H}}_{j}$
$\displaystyle\implies\textbf{A}\textbf{V}_{j}=\textbf{A}\textbf{U}(\tilde{\textbf{U}}^{*}\textbf{A}\textbf{U})^{-1}\tilde{\textbf{U}}^{*}\textbf{A}\textbf{V}_{j}+\textbf{V}_{j+1}\bar{\textbf{H}}_{j}$
$\displaystyle\implies\textbf{A}\textbf{V}_{j}=\textbf{A}\textbf{U}\textbf{B}_{j}+\textbf{V}_{j+1}\bar{\textbf{H}}_{j}$

With this extended Arnoldi relation we see that the residual update is

$\displaystyle\textbf{r}_{j}=\hat{\textbf{r}}-(\textbf{A}\textbf{U}\textbf{B}_{j}+\textbf{V}_{j+1}\bar{\textbf{H}}_{j})\textbf{y}_{j}+\textbf{A}\textbf{U}\textbf{B}_{j}\textbf{y}_{j}$
$\displaystyle=\hat{\textbf{r}}-\textbf{V}_{j+1}\bar{\textbf{H}}_{j}$

The final detail in this framework is the choice of the constraint space
$\tilde{\mathcal{V}}_{j}$ and $\tilde{\mathcal{U}}$, and as we will see this
leads to different types of augmented Krylov subspace methods. For example we
see that the choice
$\tilde{\mathcal{V}}_{j}=\textbf{A}\mathcal{K}_{j}((\textbf{I}-\Phi_{\textbf{A}\mathcal{U},\tilde{\mathcal{U}}^{\perp}})\textbf{A},\hat{\textbf{r}})$
and $\tilde{\mathcal{U}}=\textbf{A}\mathcal{U}$, in combination with imposing
$\textbf{C}=\textbf{A}\textbf{U}$ have orthonormal columns leads to the GCRO-
DR method [2]. For more details of this framework we refer the reader to [5].

##  4 Alternative formulation of the framework

The purpose of this note is to show that the projected problem (3.0.5) has an
unprojected equivalent, and from this show that by viewing the projected
problem in its unprojected form, allows the framework to easily adapt to
describe the class of so called unprojected augmented Krylov subspace methods.
We begin with the following theorem.

######  Theorem 1.

Given the projector $\Phi_{\textbf{A}\mathcal{U},\tilde{\mathcal{U}}^{\perp}}$
we have that the projected problem

Find $\textbf{t}_{j}\in\mathcal{V}_{j}$ as an approximate solution to the
problem
$\displaystyle(\textbf{I}-\Phi_{\textbf{A}\mathcal{U},\tilde{\mathcal{U}}^{\perp}})\textbf{A}\textbf{x}=(\textbf{I}-\Phi_{\textbf{A}\mathcal{U},\tilde{\mathcal{U}}^{\perp}})\textbf{r}_{0}$
subject to $\textbf{r}_{j}\perp\tilde{\mathcal{V}}_{j}$

is equivalent to the unprojected problem

Find $\textbf{t}_{j}\in\mathcal{V}_{j}$ as an approximate solution to the
problem $\displaystyle\textbf{A}\textbf{x}=\textbf{r}_{0}$ subject to
$\textbf{r}_{j}\perp(\textbf{I}-\Phi_{\textbf{A}\mathcal{U},\tilde{\mathcal{U}}^{\perp}})^{*}\tilde{\mathcal{V}}_{j}$.

###### Proof.

We will denote the residual and solution of the projected problem with the
superscript $(1)$ and that of the unprojected problem with a superscript (2).

Application of the residual constraint of the projected problem yields

$\displaystyle
0=\tilde{\textbf{V}}^{*}\textbf{r}^{(1)}_{j}=\tilde{\textbf{V}}^{*}_{j}((\textbf{I}-\Phi_{\textbf{A}\mathcal{U},\tilde{\mathcal{U}}^{\perp}})\textbf{r}_{0}-(\textbf{I}-\Phi_{\textbf{A}\mathcal{U},\tilde{\mathcal{U}}^{\perp}})\textbf{A}\textbf{V}_{j}\textbf{y}_{j}^{(1)})$
$\displaystyle\implies\textbf{y}_{j}^{(1)}=(\tilde{\textbf{V}}^{*}_{j}(\textbf{I}-\Phi_{\textbf{A}\mathcal{U},\tilde{\mathcal{U}}^{\perp}})\textbf{A}\textbf{V}_{j})^{-1}\tilde{\textbf{V}}_{j}^{*}(\textbf{I}-\Phi_{\textbf{A}\mathcal{U},\tilde{\mathcal{U}}^{\perp}})\textbf{r}_{0}$

Application of the residual constraint of the unprojected problem yields

$\displaystyle
0=\tilde{\textbf{V}}_{j}^{*}(\textbf{I}-\Phi_{\textbf{A}\mathcal{U},\tilde{\mathcal{U}}^{\perp}})\textbf{r}_{j}^{(2)}=\tilde{\textbf{V}}_{j}^{*}(\textbf{I}-\Phi_{\textbf{A}\mathcal{U},\tilde{\mathcal{U}}^{\perp}})(\textbf{r}_{0}-\textbf{A}\textbf{V}_{j}\textbf{y}_{j}^{(2)})$
$\displaystyle\implies\textbf{y}_{j}^{(2)}=(\tilde{\textbf{V}}^{*}_{j}(\textbf{I}-\Phi_{\textbf{A}\mathcal{U},\tilde{\mathcal{U}}^{\perp}})\textbf{A}\textbf{V}_{j})^{-1}\tilde{\textbf{V}}_{j}^{*}(\textbf{I}-\Phi_{\textbf{A}\mathcal{U},\tilde{\mathcal{U}}^{\perp}})\textbf{r}_{0}$

Thus $\textbf{y}^{(1)}_{j}=\textbf{y}_{j}^{(2)}$. Since both problems have the
same search space $\mathcal{V}_{j}$ we also have that in exact arithmetic
$\textbf{x}_{j}^{(1)}=\textbf{x}_{j}^{(2)}$. ∎

######  Corollary 1.1.

For the unprojected form of the problem we compute $\textbf{y}_{j}$ by solving
the following $j\times j$ linear system

$(\tilde{\textbf{V}}_{j}^{*}(\textbf{I}-\Phi_{\textbf{A}\mathcal{U},\tilde{\mathcal{U}}^{\perp}})\textbf{A}\textbf{V}_{j})\textbf{y}_{j}=\tilde{\textbf{V}}^{*}_{j}(\textbf{I}-\Phi_{\textbf{A}\mathcal{U},\tilde{\mathcal{U}}^{\perp}})\textbf{r}_{0}$
(4.0.1)

The correction for z is the same as in the original framework and is not
dependent on which form we view the problem in step one.

We now show that simply by viewing the projected problem in its unprojected
form, the framework easily adapts to describe unprojected augmented methods.
We show this using two examples. Firstly we show that the R3GMRES algorithm
arises naturally from the framework by choosing appropriate subspaces and
viewing step one of the framework as solving the unprojected problem. Using
the same procedure with correct choices of subspaces, we then derive the first
unprojected augmented FOM algorithm.

###  4.1 An unprojected augmented GMRES

The following proposition shows that by choosing subspaces corresponding to
GMRES type method (ie $\tilde{\mathcal{V}}_{j}=\textbf{A}\mathcal{V}_{j}$ and
$\tilde{\mathcal{U}}=\textbf{A}\mathcal{U}$), as well as imposing that
$\textbf{C}=\textbf{A}\textbf{U}$ have orthonormal columns, we can recover the
R3GMRES algorithm [6] . We do this by showing the linear system (4.0.1)
reduces to the same problem derived in [7, Theorem 4.1]. The result in [7] was
derived by showing that the normal equations of the R3GMRES minimization
problem can be decoupled into the equation for $\textbf{y}_{j}$ and z. The
proposition below shows this same result can be obtained by viewing the
projected problem in step one as its unprojected equivalent. We also note how
the correction for z always remains the same as in the original formulation of
the framework.

######  Theorem 2.

An unprojected augmented method which minimizes the residual over the space
$\mathcal{V}_{j}=\mathcal{K}_{j}(\textbf{A},\textbf{r}_{0})$ with the Petrov-
Galerkin constraint
$\textbf{r}_{j}\perp\textbf{A}(\mathcal{K}_{j}(\textbf{A},\textbf{r}_{0})+\mathcal{U})$,
along with imposing $\mathcal{C}=\textbf{A}\mathcal{U}$ and
$\textbf{C}^{*}\textbf{C}=\textbf{I}_{k}$ leads to the following linear system
for $\textbf{y}_{j}$.

$(\bar{\textbf{H}}_{j}^{*}\bar{\textbf{H}}_{j}-\bar{\textbf{H}}_{j}^{*}\textbf{V}^{*}_{j+1}\textbf{C}\textbf{C}^{*}\textbf{V}_{j+1}\bar{\textbf{H}}_{j})\textbf{y}_{j}=\bar{\textbf{H}}_{j}^{*}\textbf{V}_{j+1}^{*}(\textbf{I}-\textbf{C}\textbf{C}^{*})\textbf{r}_{0}$
(4.1.1)

The correction for z remains as
$\textbf{z}=\textbf{C}^{*}\textbf{r}_{0}-\textbf{B}_{j}\textbf{y}_{j}$.

###### Proof.

The orthonormality of C gives the projector
$\Phi_{\textbf{A}\mathcal{U},\tilde{\mathcal{U}}^{\perp}}=\Phi_{\mathcal{C},\mathcal{C}^{\perp}}=\textbf{C}\textbf{C}^{*}$.
Using the Arnoldi process to build a basis for
$\mathcal{K}_{j}(\textbf{A},\textbf{r}_{0})$ we obtain the usual Arnoldi
relation

$\textbf{A}\textbf{V}_{j}=\textbf{V}_{j+1}\bar{\textbf{H}}_{j}$ (4.1.2)

The linear system (4.0.1) then becomes

$\displaystyle(\tilde{\textbf{V}}_{j}^{*}(\textbf{I}-\textbf{C}\textbf{C}^{*})\textbf{A}\textbf{V}_{j})\textbf{y}_{j}=\tilde{\textbf{V}}^{*}_{j}(\textbf{I}-\textbf{C}\textbf{C}^{*})\textbf{r}_{0}$
$\displaystyle\implies((\textbf{A}\textbf{V}_{j})^{*}\textbf{A}\textbf{V}_{j}-(\textbf{A}\textbf{V}_{j})^{*}\textbf{C}\textbf{C}^{*}\textbf{A}\textbf{V}_{j})\textbf{y}_{j}=(\textbf{A}\textbf{V}_{j})^{*}(\textbf{I}-\textbf{C}\textbf{C}^{*})\textbf{r}_{0}$
$\displaystyle\implies((\textbf{V}_{j+1}\bar{\textbf{H}}_{j})^{*}(\textbf{V}_{j+1}\bar{\textbf{H}}_{j})-(\textbf{V}_{j+1}\bar{\textbf{H}}_{j})^{*}\textbf{C}\textbf{C}^{*}\textbf{V}_{j+1}\bar{\textbf{H}}_{j})\textbf{y}_{j}=(\textbf{V}_{j+1}\bar{\textbf{H}}_{j})^{*}(\textbf{I}-\textbf{C}\textbf{C}^{*})\textbf{r}_{0}$
$\displaystyle\implies(\bar{\textbf{H}}^{*}_{j}\bar{\textbf{H}}_{j}-\bar{\textbf{H}}_{j}^{*}\textbf{V}^{*}_{j+1}\textbf{C}\textbf{C}^{*}\textbf{V}_{j+1}\bar{\textbf{H}}_{j})\textbf{y}_{j}=\bar{\textbf{H}}^{*}_{j}\textbf{V}_{j+1}^{*}(\textbf{I}-\textbf{C}\textbf{C}^{*})\textbf{r}_{0}$

Thus we have derived the linear system for $\textbf{y}_{j}$ (4.1.1). It is
trivial to show that orthonormality of C also allows equation (3.0.4) to
reduce to
$\textbf{z}=\textbf{C}^{*}\textbf{r}_{0}-\textbf{B}_{j}\textbf{y}_{j}$. This
provides the update for step two. The results we derived here can be seen to
be equivalent to that shown in [7] Theorem 4.1 by setting
$\textbf{D}:=\textbf{V}^{*}_{j+1}\textbf{C}$. ∎

The full solution and the residual updates are

$\displaystyle\textbf{x}_{j}=\textbf{x}_{0}+\textbf{V}_{j}\textbf{y}_{j}+\textbf{U}\textbf{C}^{*}\textbf{r}_{0}-\textbf{U}\textbf{B}_{j}\textbf{y}_{j}$
$\displaystyle\textbf{r}_{j}=\textbf{r}_{0}-\textbf{V}_{j+1}\bar{\textbf{H}}_{j}\textbf{y}_{j}-\textbf{C}\textbf{C}^{*}\textbf{r}_{0}+\textbf{C}\textbf{B}_{j}\textbf{y}_{j}$

For full details on how to implement this version of R3GMRES, we refer the
reader to [7], wherein there is also a discussion how one can gain
implementation benefits by viewing R3GMRES in this way. An efficient MATLAB
implementation is also provided [11].

###  4.2 An unprojected augmented FOM method

We now use the framework to derive the first unprojected augmented Full
Orthogonalization Method (FOM).

######  Proposition 1.

For a FOM type method with
$\mathcal{V}_{j}=\tilde{\mathcal{V}}_{j}=\mathcal{K}_{j}(\textbf{A},\textbf{b})$
and $\tilde{\mathcal{U}}=\mathcal{U}$, with the addition of
$\textbf{C}=\textbf{A}\textbf{U}$ the unprojected method requires one to solve
the linear system

$(\textbf{H}_{j}-\textbf{V}_{j}^{*}\textbf{C}(\textbf{U}^{*}\textbf{C})^{-1}\textbf{U}^{*}\textbf{V}_{j+1}\bar{\textbf{H}}_{j})\textbf{y}_{j}=\textbf{V}^{*}_{j}(\textbf{I}-\textbf{C}(\textbf{U}^{*}\textbf{C})^{-1}\textbf{U}^{*})\textbf{r}_{0}$
(4.2.1)

The correction for z is then
$\textbf{z}=(\textbf{U}^{*}\textbf{C})^{-1}\textbf{U}^{*}\textbf{r}_{0}-\textbf{B}_{j}\textbf{y}_{j}$

###### Proof.

With these choices of search and constraint spaces we have
$\Phi_{\textbf{A}\mathcal{U},\tilde{\mathcal{U}}^{\perp}}=\textbf{C}(\textbf{U}^{*}\textbf{C})^{-1}\textbf{U}^{*}$.
Building a basis for $\mathcal{K}_{j}(\textbf{A},\textbf{r}_{0})$ using the
Arnoldi process yields the usual Arnoldi relation

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

Application of the Arnoldi relation allows one to rewrite the linear system
(4.0.1) as

$\displaystyle(\textbf{V}_{j}^{*}\textbf{V}_{j+1}\bar{\textbf{H}}_{j}-\textbf{V}_{j}^{*}\textbf{C}(\textbf{U}^{*}\textbf{C})^{-1}\textbf{U}^{*}\textbf{V}_{j+1}\bar{\textbf{H}}_{j})\textbf{y}_{j}=\textbf{V}_{j}^{*}(\textbf{I}-\textbf{C}(\textbf{U}^{*}\textbf{C})^{-1}\textbf{U}^{*})\textbf{r}_{0}$
(4.2.3)

Orthonormality of $\textbf{V}_{j}$ then yields (4.2.1). The correction for z
arises from (3.0.4) with $\tilde{\mathcal{U}}=\mathcal{U}$. ∎

The full solution and residual updates can now be written as

$\displaystyle\textbf{x}_{j}=\textbf{x}_{0}+\textbf{V}\textbf{y}_{j}-\textbf{U}\textbf{B}_{j}\textbf{y}_{j}+\textbf{U}(\textbf{U}^{*}\textbf{C})^{-1}\textbf{U}^{*}\textbf{r}_{0}$
$\displaystyle\textbf{r}_{j}=\textbf{r}_{0}-\textbf{V}_{j+1}\bar{\textbf{H}}_{j}\textbf{y}_{j}-\textbf{C}(\textbf{U}^{*}\textbf{C})^{-1}\textbf{U}^{*}\textbf{r}_{0}+\textbf{C}\textbf{B}_{j}\textbf{y}_{j}$

The method in principal allows one to recycle information from a previous
system solve and thus we refer to this method as unprojected recycled FOM or
unproj rFOM. We present an outline of the method below.

Input: $\textbf{A}\in\mathbb{C}^{n\times n}$ ,
$\textbf{x},\textbf{b}\in\mathbb{C}^{n}$ and $\textbf{U}\in\mathbb{C}^{n\times
k}$ basis for $\mathcal{U}$, and C such that
$\textbf{C}=\textbf{A}\textbf{U}$,
cycle length $j$, tolerance $\epsilon$
Result: $\textbf{x}_{j}$ such that
$\|\textbf{b}-\textbf{A}\textbf{x}_{j}\|\leq\epsilon$, updated subspace
$\mathcal{U}$

1 Algorithm
$\textbf{r}_{0}\leftarrow\textbf{b}-\textbf{A}\textbf{x}_{0}$
while  _$\|\textbf{r}_{0}\| >\epsilon$_ do

2 Build $\mathcal{K}_{j}(\textbf{A},\textbf{r}_{0})$ via the Arnoldi algorithm
to generate $\textbf{V}_{j+1}$ and $\bar{\textbf{H}}_{j}$.
Define
$\textbf{B}_{j}:=(\textbf{U}^{*}\textbf{C})^{-1}\textbf{U}^{*}\textbf{V}_{j+1}\bar{\textbf{H}}_{j}$
$\textbf{r}\leftarrow\textbf{r}_{0}-\textbf{C}(\textbf{U}^{*}\textbf{C})^{-1}\textbf{U}^{*}\textbf{r}_{0}$
Solve
$(\textbf{H}_{j}-\textbf{V}_{j}^{*}\textbf{C}(\textbf{U}^{*}\textbf{C})^{-1}\textbf{U}^{*}\textbf{V}_{j+1}\bar{\textbf{H}}_{j})\textbf{y}_{j}=\textbf{V}^{*}_{j}\textbf{r}$
$\textbf{x}\leftarrow\textbf{x}+\textbf{V}_{j}\textbf{y}_{j}-\textbf{U}\textbf{B}_{j}\textbf{y}_{j}+\textbf{U}(\textbf{U}^{*}\textbf{C})^{-1}\textbf{U}^{*}\textbf{r}_{0}$
$\textbf{r}_{0}\leftarrow\textbf{r}-\textbf{V}_{j+1}\bar{\textbf{H}}_{j}\textbf{y}_{j}+\textbf{C}\textbf{B}_{j}\textbf{y}_{j}$
Update $\mathcal{U}$ and $\mathcal{C}$ such that
$\mathcal{C}=\textbf{A}\mathcal{U}$.

3 end while

4

Algorithm 3 unprojected Recycled FOM (unprojected rFOM($m$))

### Numerical Experiment

We have implemented the unproj rFOM algorithm in MATLAB to demonstrate its
effectiveness. Below we demonstrate results of applying the algorithm to a
Wilson Dirac matrix arising in Lattice QCD. The augmentation subspace
$\mathcal{U}$ was chosen to be the eigenvectors corresponding to the $k$
smallest eigenvalues obtained from a closely related matrix.

Figure 1: An example run of the unprojected recycled FOM algorithm (unproj
rFOM) compared to the standard FOM algorithm of Saad [8]. This test was run on
Wilson Dirac matrix of size $49,152\times 49,152$ arising in Lattice QCD. The
matrix corresponds to a $8^{4}$ lattice. The cycle length was $m=40$ and
different values of $k$ were chosen to demonstrate how larger $k$ allows for
faster convergence (of course, with additional computational expense).

##  5 Conclusion

In this paper we have taken an alternative view of the general framework
describing augmented Krylov subspace methods described in [5]. We have shown
that by viewing the projected problem arising in this framework in terms of
its unprojected equivalent, the framework easily and naturally extends to
describe unprojected augmented Krylov subspace methods. We have demonstrated
that for a GMRES type method, the framework reduces to the R3GMRES algorithm
[6]. We have further developed an unprojected augmented FOM algorithm and have
demonstrated its effectiveness as a recycle method.

There is an active area of interest in unprojected augmented Krylov subspace
methods, particularly for solving ill posed inverse problems [7, 6], and
viewing the original framework allows one to easily derive new such methods as
we have demonstrated in this paper.

## References

  * [1] André Gaul, Martin H Gutknecht, Jorg Liesen and Reinhard Nabben  “A framework for deflated and augmented Krylov subspace methods”  In _SIAM Journal on Matrix Analysis and Applications_ 34.2 SIAM, 2013, pp. 495–518
  * [2] Michael L Parks et al.  “Recycling Krylov subspaces for sequences of linear systems”  In _SIAM Journal on Scientific Computing_ 28.5 SIAM, 2006, pp. 1651–1674
  * [3] Martin H Gutknecht  “Deflated and augmented Krylov subspace methods: A framework for deflated BiCG and related solvers”  In _SIAM Journal on Matrix Analysis and Applications_ 35.4 SIAM, 2014, pp. 1444–1466
  * [4] André Gaul  “Recycling Krylov subspace methods for sequences of linear systems”, 2014
  * [5] Kirk M. Soodhalter, Eric Sturler and Misha Elena Kilmer  “A survey of subspace recycling iterative methods”  In _GAMM‐Mitteilungen_ 43, 2020
  * [6] Yiqiu Dong, Henrik Garde and Per Christian Hansen  “R3GMRES: including prior information in GMRES-type methods for discrete inverse problems”  In _Electronic Transactions on Numerical Analysis_ 42 Kent State University Institute of Computational Mathematics, 2014, pp. 136–146
  * [7] Kirk M Soodhalter  “A note on augmented unprojected Krylov subspace methods”  In _arXiv preprint arXiv:2106.10050_ , 2021
  * [8] Yousef Saad  “Krylov subspace methods for solving large unsymmetric linear systems”  In _Mathematics of computation_ 37.155, 1981, pp. 105–126
  * [9] Yousef Saad and Martin H Schultz  “GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems”  In _SIAM Journal on scientific and statistical computing_ 7.3 SIAM, 1986, pp. 856–869
  * [10] Yousef Saad  “Iterative methods for sparse linear systems”  SIAM, 2003
  * [11] Kirk M Soodhalter  “kirkmsoodhalter/r3gmres-simplified: R3gmres simplified implementation”  URL: https://doi.org/10.5281/ZENODO.4975990