# Improving Antibody Humanness Prediction using Patent Data

Talip Uçar  Aubin Ramon  Dino Oglic  Rebecca Croasdale-Wood  Tom Diethe
Pietro Sormanni

###### Abstract

We investigate the potential of patent data for improving the antibody
humanness prediction using a multi-stage, multi-loss training process.
Humanness serves as a proxy for the immunogenic response to antibody
therapeutics, one of the major causes of attrition in drug discovery and a
challenging obstacle for their use in clinical settings. We pose the initial
learning stage as a weakly-supervised contrastive-learning problem, where each
antibody sequence is associated with possibly multiple identifiers of function
and the objective is to learn an encoder that groups them according to their
patented properties. We then freeze a part of the contrastive encoder and
continue training it on the patent data using the cross-entropy loss to
predict the humanness score of a given antibody sequence. We illustrate the
utility of the patent data and our approach by performing inference on three
different immunogenicity datasets, unseen during training. Our empirical
results demonstrate that the learned model consistently outperforms the
alternative baselines and establishes new state-of-the-art on five out of six
inference tasks, irrespective of the used metric.

Machine Learning, ICML

##  1 Introduction

Proteins are naturally occurring large molecules, composed of one or more
chains of amino acids. They are present in all living organisms and include
many essential biological compounds such as enzymes and antibodies. An
antibody is a type of protein that binds specifically to proteins called
antigens, typically found on the surface of pathogens such as bacteria,
viruses, and other foreign substances, including cancer cells. Immune system
protects the body from them by having a large number of antibodies in
circulation throughout the body, capable of binding specifically to different
antigens. Monoclonal antibodies are manufactured in laboratories as clones or
copies of a single antibody. They are a significant class of
biopharmaceuticals, often used for therapeutic purposes, and have contributed
significantly to the overall sales of biopharmaceutical products (Ecker et
al., 2015). Therapeutics based on monoclonal antibodies are typically produced
from non-human sources and can lead to an undesired immune response, referred
as _immunogenicity_ , which is one of the major causes of attrition in the
development of monoclonal antibody therapeutics (Vandivort et al., 2020).
Thus, many therapeutic monoclonal antibodies often go through a _humanization_
process, which involves modifying the antibody to make it more similar to
human antibodies, reducing the likelihood of triggering an immune response.
This process must be able to identify the humanness of an antibody with little
to no errors as there are significant costs associated with attrition in the
later stages of the drug development process (Marks et al., 2021).

Motivated by the importance of antibody humanness prediction, we present a
Self-supervised learning framework designed to extract informative and
contextually rich representations of amino acids and antibodies from Patented
Antibody Database (SelfPAD). The Patented Antibody Database (PAD) is a rich
and diverse source of information with functional characterization of
antibodies, consisting of approximately $300,000$ unpaired sequences across
seven different species (Krawczyk et al., 2021). Our approach to leverage the
PAD for representation learning and immunogencity prediction departs from
conventional approaches, pre-dominantly focused on leveraging information from
Observed Antibody Space (OAS). The latter is the database with more than
billion different antibody sequences sourced from more than $80$ different
studies (Olsen et al., 2022a). What lacks in the OAS is functional
characterization of antibodies which limits representation learning to self-
supervised approaches focused on sequence completion tasks (Ruffolo et al.,
2021; Olsen et al., 2022b; Shuai et al., 2023). In contrast, the PAD offers
noisy association of sequences with function via more than $16,000$ different
patent filings that can be translated into more than $1,000$ tasks, one for
each set of potential targets (please see Section 4.1). This allows for multi-
stage, multi-loss supervised learning, where a model is optimized with two
different loss functions during the process. More specifically, we devote the
initial stage to self-supervised representation learning to learn the
representations of amino acids. We then switch the loss function to focus the
final epochs entirely on a predictive task of interest, e.g. humanness
prediction.

We pose the initial stage focused on representation learning from patents as a
weakly-supervised contrastive-learning problem, where each antibody sequence
is associated with possibly multiple noisy identifiers of function (i.e.
potential targets) and the objective is to learn a latent representation that
groups them according to their patented objectives. In the first instance,
this allows for learning of an antibody representation directly by the
contrastive encoder that is informed pre-dominantly by the specific binding to
different targets. One of our underlying assumptions is that the amino acid
embeddings acquired through our proposed framework carry information about the
functional and developability traits of antibodies as the patent filing is
typically representative of the drug discovery process where the initial goal
is to identify a potent binder, followed by subsequent improvements focused on
refining its developability properties, encompassing factors such as
immunogenicity, aggregation, viscosity etc. To demonstrate the utility of
patent data for improving the immunogenicity prediction, we continue training
the contrastive encoder to predict the humanness score of a given antibody
sequence. We explore two different avenues for achieving transfer learning
between the stages, the first via the pre-trained encoder itself and the
second one via the learned embeddings of individual amino acids. During this
stage of the training process, we also add a multi-layer perceptron (MLP) on
top of the contrastive encoder and continue training on the patent data using
the cross-entropy loss function. We evaluate the effectiveness of the learned
immunogenicity model by performing inference on three different immunogenicity
datasets published as part of prior work (see Sections 2 and 4.1). Our
empirical analysis demonstrates that the final model learned from patent data
achieves the new state-of-the-art performance on five out of six inference
tasks, irrespective of the used metric. SelfPAD is available at:
https://github.com/AstraZeneca/SelfPAD

##  2 Related Work

(a) Pre-training

(b) Fine-tuning

Figure 1: SelfPAD Framework: a) Pre-training: The training process that
considers sequences associated with the same potential targets as positive
pairs. b) Fine-tuning: The architecture of the transformer-based sequence
encoder used in pre-training as well as MLP, swap noise and shift noise used
only during fine-tuning.

Various _humanness_ scores have been developed to measure the human-like
characteristics of humanized antibody sequences and identify potential risks
for immunogenicity. Among these, two of the most prominent methods determine
humanness by considering sequence identity with a reference library of human
sequences, either averaged across all sequences (Z-score, Abhinandan & Martin,
2007) or across the closest $20$ sequences (T20 score, Gao et al., 2013). Both
of these approaches provide a single score for the entire sequence and rely on
small reference libraries, restricting the diversity of designed therapeutics.
Recently, novel methods incorporating machine learning (ML) have been proposed
that leverage patterns associated with human and non-human sequences to
predict the humanness of an antibody. The most popular approaches along these
lines are Hu-mAb (Marks et al., 2021), OASis (Prihoda et al., 2022), and
AbNatiV (Ramon et al., 2024).

OASis is an approach quantifying the humanness of an antibody sequence based
on peptide $k$-mers. It works by first segmenting a sequence into overlapping
$9$-mer peptides and then querying these fragments against the Observed
Antibody Space (OAS) (Olsen et al., 2022a) to estimate their prevalence within
the human population. To comprehensively evaluate and compare humanness at the
entire antibody level, it introduces the _OASis identity score_ , which is
calculated at the antibody level as the fraction of its $k$-mer peptides that
surpass a user-defined prevalence threshold, where the threshold dictates the
required fraction of the human population for a given peptide to be deemed
human. OASis offers four predefined thresholds, capturing different stringency
levels: _i_) loose ($\geq 1\%$ subjects), _ii_) relaxed ($\geq 10\%$
subjects), _iii_) medium ($\geq 50\%$ subjects), and _iv_) strict ($\geq 90\%$
subjects). For instance, at the relaxed threshold, a peptide is classified as
human if found in at least $10\%$ of subjects. Accordingly, the OASis relaxed
identity score for an antibody is computed as the fraction of its peptides
found in at least $10\%$ of subjects. It is important to note that OASis may
not perform effectively for antibodies with long CDR3 loops (Prihoda et al.,
2022) since CDR3 loops are highly varied regions.

AbNatiV is based on a vector-quantized variational auto-encoder (VQ-VAE) (Van
Den Oord et al., 2017) trained with masked unsupervised learning on aligned
variable domain (Fv) sequences from curated native immune repertoires. Four
distinct models are trained on Fv sequences corresponding to _i_ -_iii_) human
heavy (VH), kappa (VKappa) and lambda (VLamda) chains, and _iv_) camelid heavy
chain single-domains (VHH or nanobodies). It is specifically designed for
quantifying the humanness of antibodies and nanobodies. Trained on curated
datasets of human antibodies or camelid nanobodies, it predicts the likelihood
of a given sequence belonging to the distribution of immune-system-derived
antibodies. It can be used for both designing antibodies and nanobodies, and
predicting the likelihood of immunogenicity. One potential shortcoming of
AbNatiV is that it relies on four separate models, one for each chain type,
and, thus, fails to account for the entirety of information required for
quantifying the humanness property.

Hu-mAb is another approach that utilizes the OAS database and trains multiple
random forest (RF) classifiers to distinguish between human and non-human
sequences, each tailored for a specific human V gene type. It can be used to
infer the _humanness_ of sequences as well as for humanization of sequences.
Similar to AbNativ, relying on multiple models is one of its potential
shortcomings. It is also known to perform worse with light chains (Marks et
al., 2021).

When it comes to machine learning aspects of our work, it is related to setups
involving self-supervised pre-training and task-specific fine-tuning. The
paradigm of pre-training, whether supervised or unsupervised, and subsequent
fine-tuning has proven to be a potent approach with promising results across
various domains such as natural language processing (Mikolov et al., 2013;
Pennington et al., 2014; Devlin et al., 2018; Liu et al., 2019; Conneau et
al., 2019; Joulin et al., 2016; Collobert & Weston, 2008), computer vision and
audio domains (Chen et al., 2020a; Grill et al., 2020; Oord et al., 2018;
Caron et al., 2020; He et al., 2020; Falcon & Cho, 2020; Chen et al., 2020b).
This progress was mainly enabled by taking advantage of spatial, semantic, or
temporal structure in the data through data augmentation (Chen et al., 2020a),
pretext task generation (Devlin et al., 2018) and using inductive biases
through architectural choices (e.g., convolutional layers for images). In this
particular instance, we leverage contrastive learning (Chen et al., 2020a)
tailored for the multi-task setting at the initial epochs of the training
process, and this stage can be considered a pre-training step. The second
phase of the training process can be viewed as a form of fine-tuning. Unlike
the conventional approach found in the literature, this fine-tuning occurs on
the original dataset with the objective of learning to quantify the humanness
of an antibody sequence. We should note that a well-curated patent dataset
could facilitate the training of a multi-label, multi-task classifier.
However, the identifiers of function currently available are noisy and may not
be highly reliable. This is often due to a tendency toward risk aversion in
safeguarding intellectual property, potentially resulting in patent filings
containing sequences with varying functions. An alternative strategy for pre-
training might involve disregarding the function and concentrating on self-
supervised learning, which typically involves sequence completion tasks.
However, we leave this as a potential avenue for future exploration. Finally,
it’s important to emphasize that patent data serves as a unique information
source for assessing the humanness of antibodies. Existing works, including
those discussed here, predominantly rely on the OAS database. However, it is
essential to note that the OAS database exhibits bias and lacks the necessary
diversity to comprehensively characterize the entire population of human
antibodies.

##  3 Method

In this section, we describe our approach for learning an effective model for
predicting the humanness of an antibody sequence. As described in Section 1,
we perform a multi-stage training process with two different loss functions.
There are two main steps characterizing our approach: _i_) learning an
antibody sequence encoder by pre-training on patent data using weakly
supervised constrastive learning, and _ii_) learning an antibody property
prediction model by fine-tuning the encoder on the same dataset using the
cross-entropy loss function and humanness as the task.

Pre-training The encoder consists of two components: _i_) a core transformer-
based encoder that performs one-to-one mapping between amino acids from the
antibody sequence and their context-dependent latent space embeddings and
_ii_) a projection layer that maps the latent space representation of an
antibody sequence to a fixed dimensional vector, independent of the sequence
length. The transformer block has four layers and each one is implemented such
that the layer normalization is applied prior to the self-attention and
feedforward blocks (e.g., see Radford et al., 2019). We use a maximum sequence
length of $180$ and the dimension of $32$ for the embedding and hidden layers.
At the input, we first encode each amino acid in the sequence by using two
trainable embedding layers, based on amino-acid type and position in the
sequence, i.e., its positional encoding. We sum these two embedding vectors
and concatenate the resulting embedding with additional $18$ features that are
pre-computed based on the biophysical properties of amino acids such as their
hydrophobicity, charge, volume etc. (please see Table 7 in the appendix for
the full list of properties). We also inject a $\sigma$-scaled Gaussian noise
(with $\sigma=0.01$) to the embeddings to learn a robust amino acid
representation during training. We do not use any alignment for the sequences
and do not explicitly condition the model with respect to the type of chain or
the origin of species, i.e., we use each sequence as given. Since we use
unigrams, the vocabulary size is twenty three to account for twenty amino
acids as well as the tokens used for beginning and end of the sequence, and
padding (#, $, and * respectively). The projection layer includes a linear
layer, followed by GELU activation and a final linear layer, where both linear
layers have a dimension of $128$. We pre-trained the model with a batch size
of $100$ for $1000$ epochs (see Figure 5 in the appendix for details).

Figure 1a depicts our contrastive learning setup. Antibody sequences are
grouped by their association with potential function identifiers as specified
in patents. Let $K$ be the total number of targets available. We uniformly
sample $k$ targets without replacement from this set in every iteration during
training, i.e.,

$\displaystyle T=\\{t_{1},...,t_{k}\\},\quad t_{i}\sim\mathcal{U}(1,K),\quad
i=1,..,k,$ (1)

where $t_{i}$ denotes the index of the $i$-th sampled target. This process
ensures that each element in the set of $K$ targets has an equal probability
of being selected, resulting in the creation of a representative subset $T$ of
size $k$. Then, we uniformly sample two sequences (positive pairs) for each of
$k$ targets, resulting in a minibatch size of $2k$:

$\displaystyle S=\\{(s_{i1},s_{i2})\mid
s_{i1},s_{i2}\sim\mathcal{U}(G_{t_{i}}),i=1,..,k\\},$ (2)

where $(s_{i1},s_{i2})$ represents a pair of sequences from a group of
sequences $G_{t_{i}}$ associated with the target $t_{i}$. This approach
ensures that for each target $t_{i}$, two sequences are uniformly chosen,
thereby generating a set of representative samples. Given a positive pair of
sequences $(s_{i1},s_{i2})$, we treat the remaining $2(k-1)$ sequences within
a minibatch as negative examples. As shown in Figure 1a, we obtain embeddings
$(z_{i1},z_{i2})$ for each pair of sequences by using a shared encoder and a
projection network. Then, we compute the normalized temperature-scaled cross
entropy loss (Chen et al., 2020a) for a positive pair
$(\bm{z}_{i1},\bm{z}_{i2})$ as

$\displaystyle\begin{aligned}
&\mathcal{L}_{c}=\frac{1}{2k}\sum_{i=1}^{k}\left[l(\bm{z}_{i1},\bm{z}_{i2})+l(\bm{z}_{i2},\bm{z}_{i1})\right]\quad\text{with}&\\\
&l(\bm{z}_{i1},\bm{z}_{i2})=-\log\frac{\exp(\psi(\bm{z}_{i1},\bm{z}_{i2})/\tau)}{\sum_{j=1,j\neq
i}^{2k}\exp(\psi(\bm{z}_{i1},\bm{z}_{j2})/\tau)}&\end{aligned}$ (3)

where $\psi(.)$ is the dot product similarity score, $\mathcal{L}_{c}$ is
total contrastive loss, and $l(\bm{z}_{i1},\bm{z}_{i2})$ is the loss function
for a positive pair $(\bm{z}_{i1},\bm{z}_{i2})$.

(a) Species

(b) Potential targets

(c) Clusters of chain types

Figure 2: After pre-trainnig: t-SNE visualisation of antibodies in the test
set based on their a) Origin and b) Potential targets (showing only small
subset of targets for better readibility). After fine-tuning the model to
predict for humanness: c) t-SNE plot of heavy, kappa, lambda and vhh chains
from different species. H-Heavy=Heavy chain from humans, O-Heavy: Heavy chain
from other species.

Fine-tuning We add a multi-layer perceptron (MLP) on top of the transformer-
based encoder learned in the pre-training step. MLP consists of blocks with a
residual connection, where each block has a linear layer with dimension
$128\times 512$ followed by GELU activation, batch normalisation, dropout
(with $p=0.3$) and another linear layer with dimension $512\times 128$. The
final linear layer is paired with the softmax activation to give a binary
classifer. To learn a robust representation of amino acids and antibodies, we
also apply two types of noise during the training: _i_) shift noise: half of
the time, we shift the starting position of a sequence by a randomly sampled
number from the range $[0,10]$, and mask tokens are added to the beginning,
shifting the sequence to the right. _ii_) swap noise: before inputting the
resulting embedding into the MLP, we incorporate swap noise (Ucar et al.,
2021), which involves randomly replacing selected embedding entries with
values sampled from the same column.

We explored different strategies for the fine-tuning phase of the training
process, where for example, one re-initializes the final four layers of
encoder (effectively keeping only the learned amino acid embeddings and
position encoding fixed) and fine-tunes by applying different learning rates
to each layer of the contrastive encoder as well as the MLP blocks. More
specifically, we increase the learning rate from low $0.0001$ to high $0.001$
as we move from lower to upper layers in the encoder. Alternatively, we also
investigated the performance of the model that continues training of the
contrastive encoder (i.e., no re-initialization of layers) on the humanness
task during the fine-tuning epochs (including the extra MLP block). In the
second phase of the training, we used cross-entropy loss with label smoothing
(configured to $0.5$) and a batch size of $512$ for $25$ epochs.

##  4 Experiments

In this section, we present our empirical results and provide evidence for the
effectiveness of the approach in improving the humanness prediction of
antibodies. We start with a high-level description of the used datasets and
data-processing performed prior to training. Then, we summarize our empirical
observations for the pre-training and fine-tuning steps of the multi-stage
training process and illustrate the effectiveness of the approach relative to
baselines.

###  4.1 Data Sources and Pre-processing

Patent data for antibodies can be a potentially valuable source of information
for learning useful representations from their sequences. The primary purpose
of patents is to provide legal protection rather than to serve a scientific
knowledge, hence such database can be noisy. For example, although patented
antibodies are tagged with potential targets, their relationship is not always
clear. However, despite of its low signal-to-noise ratio, it can contain
valuable information on relationships between family of antibody sequences and
their potential targets. To illustrate its utility, we use Patented Antibody
Database (PAD) (Krawczyk et al., 2021) that contains $16,526$ patent families
from major jurisdictions such as US Patent and Trademark Office and World
Intellectual Property Organization. These families hold $\sim 290k$ unique
antibody chains (unpaired heavy and light chains as well as nanobodies) that
are compiled in the PAD. Sequences in the database are manually tagged with
potential target(s). If a sequence is associated with more than one target
(e.g. ’il13’, ’il4’), we use the set of targets as the label (e.g. [’il13’,
’il4’]). We group sequences based on their association with potential targets.
Sequences in each group constitute our positive samples for the corresponding
target. We then randomly split the sequences into training ($\sim 260k$) and
test sets ($\sim 29k$) by using stratified sampling to keep the proportion of
antibodies from humans and non-humans same in both sets. Since the PAD also
includes sequences from other publicly available datasets, we excluded $553$
therapeutics and $217$ immunogenicity data (Prihoda et al., 2022), along with
$25$ humanization data (Marks et al., 2021), from the training set to avoid
data leakage, considering that we utilize these datasets for evaluation
purposes. Additionally, we excluded all sequences exhibiting at least $95\%$
sequence similarity with any sequence in the latter three datasets. The
minimum Levenshtein distance of sequences in the three evaluation datasets
from those in the training set of the PAD is presented in Table 1. Moreover,
during pre-training, we drop any sequence that is not tagged with any target
as well as any target that has less than 20 sequences in its group. After
filtering, we end up with $1063$ potential targets and $\sim 136k$ sequences
in the training set, which we then use to pre-train the encoder. In fine-
tuning the pre-trained encoder for humanness prediction, we utilize all
available training set sequences. The training set is then split into $90\%$
for training and $10\%$ for validation, and the model is fine-tuned with
cross-entropy loss for $25$ epochs.

Table 1: Minimum Levenshtein distance between evaluation datasets and the PAD
training set. Abbreviations are; LD: Levenshtein Distance, 553 T-H/T-L: 553
Therapeutics data Heavy/Light chains, 217 Im.-H/Im.-L: 217 Immunogenicity data
Heavy/Light chains, 25 Hum.-P/Hum.-H 25 Parental/Humanized sequences from
humanization data.

| 553 T-H | 553 T-L | 217 Im.-H | 217 Im.-L | 25 Hum.-P | 25 Hum.-H
---|---|---|---|---|---|---
LD | 8 | 6 | 4 | 8 | 6 | 9

(a) PAD test set

(b) 553 therapeutics

(c) 217 ADA data

(d) 25 humanization data

Figure 3: Human versus Nonhuman sequences: PCA of embeddings from FT-SelfPAD
for a) The test set of PAD, b) 553 therapeutics data from OASis, c) 217
immunogenicity data from OASis, d) 25 humanization data from Hu-mAb, where 25
therapeutic antibodies with immunogenicity problem (i.e. parental sequence) as
well as their humanized versions, plotted together with rest of the sequences
in the test set from the PAD (shown in blue).

553 Therapeutics is a dataset from Prihoda et al. (2022), also used in OASis,
and originally sourced from IMGT mAb DB (Poiron et al., 2010). After removing
four antibodies (letolizumab, lulizumab pegol, placulumab and glofitamab), we
are left with $449$ paired sequences, consisting of $197$ human, $226$
humanized, $63$ chimeric, $41$ humanized/chimeric, $13$ murine, $6$ caninized,
and $3$ felinized sequences. For humanness prediction evaluation, we classify
the positive class as containing human and humanized sequences ($423$
sequences), while the rest are labeled negative.

217 immunogenicity refers to the dataset obtained from Prihoda et al. (2022)
and also used in OASis. It contains reported ADA responses from $217$ paired
sequences and is originally curated in Marks et al. (2021). We evaluate our
model by scoring $217$ therapeutics for humanness. Sequences with
immunogenicity $\geq 10$ are classified as nonhuman while the ones with $<10$
are as human.

25 humanization data refers to the dataset with $7$ pairs of parental-
humanized sequences that were first curated by Clavero-Álvarez et al. (2018)
and expanded to $25$ pairs by Marks et al. (2021). It has also been used in
the evaluation of Hu-mAb. We use these $25$ experimentally validated pairs to
predict reported ADA responses in the form of humanness prediction. We
classify parental sequences with immunogenicity issues as non-human and their
corresponding humanized versions as human.

###  4.2 Empirical Analysis of the Pre-training Stage

Figure 2 shows the clustering of sequences from test set of the PAD, based on
their origin and potential targets (Figure 2a-b). We see that the sequences
from same species as well as those associated with same potential targets tend
to cluster together. Note that when training our models, sequence embeddings
are generated by aggregating over the amino acids to prevent the model from
memorizing the properties that depend on the sequence length (though it might
be useful to learn this dependency in some context).

Table 2: The PAD test set: Summary of performance for humanness prediction.
Please note that the sequences from alpaca are excluded when computing metrics
for AbNatiV and Hu-mAb.

Model | F1 | Recall | Precision | Accuracy | ROC AUC | PR AUC
---|---|---|---|---|---|---
FT-SelfPAD | 97.59 | 97.36 | 97.83 | 96.58 | 99.12 | 98.54
AbNativ | 91.59 | 90.68 | 92.52 | 87.54 | 91.52 | 96.49
OASis | 88.67 | 86.45 | 91.00 | 83.45 | 88.30 | 95.57
Hu-mAb | 70.35 | 55.22 | 96.87 | 64.72 | 74.83 | 93.01

###  4.3 Empirical Analysis of the Fine-tuning Stage

A humanness score should be able to distinguish between human antibodies and
antibodies from other species. In particular, it should enable doing so for
therapeutic antibodies since those are the primary subjects of humanness
analysis. The goal of humanness evaluation is to capture and reduce the
immunogenicity risk. If antibodies are not structurally similar to natural
human antibodies, there is an increased risk of the immune system recognizing
them as foreign and producing antibodies against therapeutic antibodies i.e.
”anti-drug antibodies” (ADAs) or ”immunogenicity”. Therefore, we fine-tune and
evaluate SelfPAD to predict the humanness of sequences. We refer to the fine-
tuned model as FT-SelfPAD. We should note that, in the machine learning
literature, the term ”fine-tuning” is used to mean supervised training of a
model, which is pre-trained on a large corpus of data (typically in
unsupervised framework), on a small dataset. In our context, fine-tuning
corresponds to further training of pre-trained model on the same dataset, but
using a supervised loss instead of a contrastive loss so that the model can
learn to distinguish between human and nonhuman sequences. Once trained, we
benchmark FT-SelfPAD against three strong baselines, namely Hu-mAb (Marks et
al., 2021), OASis (Prihoda et al., 2022) and AbNatiV (Ramon et al., 2024), by
using four datasets: i) 29k unpaired sequences from the test set of the PAD,
ii / iii) paired sequences from 553 therapeutics and 217 immunogenicity data
(Prihoda et al., 2022) and iv) paired sequences from 25 humanization data
(Marks et al., 2021).

29k sequences from the test set of PAD Once the model is fine-tuned for
humanness prediction, we extract embeddings of the test set, and cluster
sequences based on their chain type and origin. As shown in Figure 2c, the
sequences from different origins and chain types make up distinct clusters.
For examples, heavy chains from humans (orange in Figure 2c) cluster together,
but are separate from the heavy chains of other species (purple). Similarly,
nanobodies (i.e. VHH sequences) constitute a distinct cluster away from the
rest of the sequences. This result indicates that our model should be able to
distinguish the origin and chain type of a sequence. Moreover, we note that
lambda chains from both humans and other species cluster together, hinting
that our model might have more difficulty in distinguishing the origin of
lambda chains. We also plot sequences based on their origin by using PCA of
embeddings in Figure 3a, in which we clearly see that FT-SelfPAD can
distinguish the sequences from humans, alpaca and other species. Table 2
summarizes the performance of FT-SelfPAD and other baselines for humanness
prediction and shows that FT-SelfPAD performs the best across all metrics. We
should note that AbNatiV (Ramon et al., 2024) has four models, one for each
chain type of heavy, lambda, kappa and VHH. Thus, we used AbNatiV’s
corresponding models to predict the humanness (with a threshold of 0.8) of
heavy, kappa and lambda chains from human and other species in the test set
after removing VHH sequences. For OASis, we evaluate and compare humanness by
using the OASis identity score and prevalence threshold of ’relaxed’. We
evaluated OASis at three different OASis identity threshold: 0.5, 0.7 and 0.8.
We obtained the best result by using identity threshold of 0.7, which we also
use in the rest of the experiments. Finally, Hu-mAb relies on separate
classifiers, each of which is trained for each human V gene type. It also
computes humanness score based on each chain type i.e. heavy or light chain.
If it identifies any V gene in the sequence as human, we consider the entire
sequence as human. Hence, we scored heavy and light chains from humans and
other species after removing VHH sequences. In our experiments, we noticed
that Hu-mAb performs poorly on light chains from humans, an observation that
is also reported in (Marks et al., 2021), resulting in a poor performance
overall.

Table 3: 553 therapeutics: Summary of performance for paired sequences (top) &
single chains (bottom). Please note that we consider both human and humanized
sequences as human.

Model | F1 | Recall | Precision | Accuracy | ROC AUC | PR AUC
---|---|---|---|---|---|---
Paired | FT-SelfPAD | 93.33 | 99.29 | 88.05 | 89.07 | 94.47 | 98.16
AbNativ | 91.03 | 95.98 | 86.57 | 85.43 | 90.02 | 96.73
OASis | 91.97 | 93.91 | 90.11 | 87.34 | 89.48 | 96.31
Hu-mAb | 91.89 | 95.04 | 88.94 | 87.07 | 77.68 | 93.90
Single | FT-SelfPAD | 92.00 | 93.74 | 90.32 | 87.43 | 91.84 | 97.02
AbNativ | 87.93 | 87.00 | 88.89 | 81.60 | 86.91 | 95.48
OASis | 88.14 | 87.00 | 89.32 | 81.97 | 85.42 | 94.88
Hu-mAb | 62.05 | 47.64 | 88.96 | 55.10 | 63.90 | 88.47

Table 4: 217 immunogenicity data: Summary of performance for paired sequences
(top) & single chains (bottom). Similar to (Prihoda et al., 2022), we consider
any sequence with immunogenicity $\geq 10$ as nonhuman and with $<10$ as
human.

Model | F1 | Recall | Precision | Accuracy | ROC AUC | PR AUC
---|---|---|---|---|---|---
Paired | FT-SelfPAD | 87.50 | 89.63 | 85.47 | 80.65 | 82.86 | 93.79
AbNativ | 87.91 | 90.85 | 85.14 | 81.11 | 79.90 | 92.09
OASis | 85.20 | 85.98 | 84.43 | 77.42 | 78.78 | 92.32
Hu-mAb | 88.96 | 90.85 | 87.13 | 82.95 | 74.67 | 92.45
Single | FT-SelfPAD | 83.75 | 81.71 | 85.90 | 76.04 | 78.94 | 90.44
AbNativ | 83.12 | 79.57 | 77.86 | 75.58 | 77.86 | 91.25
OASis | 82.94 | 80.79 | 85.21 | 74.88 | 76.86 | 90.75
Hu-mAb | 59.72 | 45.43 | 87.13 | 53.69 | 62.34 | 86.90

553 therapeutics data We process the dataset such that 197 human and 226
humanized sequences are labeled as human while the rest is labelled as
nonhuman with the assumption that humanized sequences is more likely to be
similar to human sequences. We evaluate the models under two scenarios: i)
Paired: We measure the humanness prediction of a paired sequence. If at least
one of two chains (light and heavy) in a paired sequence is predicted as
human, we consider the paired sequence as human. ii) Single chain: We measure
the performance of models for each chain. Table 3 summarizes the performance
of all models, showing that FT-SelfPAD performs the best for both scenarios.

217 immunogenicity data It contains reported ADA responses from 217 sequences.
Similar to (Prihoda et al., 2022), we label sequences with immunogenicity
$\geq 10$ as nonhuman while the ones with $<10$ as human. However, we should
emphasize that there is a difference between humanness and immunogenicity, and
our evaluation of humanness is only a proxy in this context. For example, some
fully human antibodies given at high dosage for extended periods of time can
end up being more immunogenic than humanised antibodies administered once at
low dose. Keeping this in mind, we list our results in Table 4, showing that
FT-SelfPAD gives the best performance for the most important metrics, namely
ROC AUC and Precision-Recall (PR) AUC, for paired sequences while it performs
best across all metrics, except PR-AUC, for the single chains.

25 humanization data The data consists of 25 parental sequences with
immunogenicity issues and their corresponding humanized versions. We label
parental sequences as nonhuman while humanized ones as human. As shown in
Table 5, FT-SelfPAD classifies all sequences correctly. Moreover, we plot both
parental sequences and humanized counterparts by using a PCA on the their
embeddings from FT-SelfPAD as shown in Figure 3d. We also overlay the
sequences from the test set of PAD to give some context. By comparing to the
Figure 3a, we see that parental sequences show up in clusters of other species
while the humanized sequences appear among the sequences from human origin.
Hence, our model is correctly able to identify the sequences with
immunogenicity issues. We made similar comparisons for both 553 therapeutics
and 217 immunogenicity data in Figure 3b-c, resulting in a similar
observation.

Table 5: 25 humanization data: Summary of performance for 25 sequences with
immunogenicity problem and their corresponding humanized versions (Marks et
al., 2021) for both paired sequences (top) & single chains (bottom).

Model | F1 | Recall | Precision | Accuracy | ROC AUC | PR AUC
---|---|---|---|---|---|---
Paired | FT-SelfPAD | 100 | 100 | 100 | 100 | 100 | 100
AbNativ | 92.00 | 92.00 | 92.00 | 92.00 | 99.52 | 99.53
OASis | 96.00 | 96.00 | 96.00 | 96.00 | 98.00 | 97.52
Hu-mAb | 96.15 | 100 | 92.59 | 96.00 | 96.00 | 96.30
Single | FT-SelfPAD | 100 | 100 | 100 | 100 | 100 | 100
AbNativ | 85.71 | 78.00 | 95.12 | 87.00 | 96.12 | 95.98
OASis | 81.72 | 76.00 | 88.37 | 83.00 | 93.14 | 92.36
Hu-mAb | 64.94 | 50.00 | 92.59 | 73.00 | 73.00 | 83.80

Interpretability To enhance the model’s utility in antibody engineering,
humanness evaluation should offer interpretable scores by highlighting the
impact of individual residues or sequence segments. Improved interpretability
includes providing explanations of humanness scores through comparisons with
closely related reference sequences. We tested FT-SelfPAD’s ability to infer
positions for humanizing sequences by analyzing attention in 25 humanization
data sequences. To identify potential candidates for humanization, we
subtracted attention scores (extracted from the last layer of encoder) of
humanized sequences from their corresponding parental sequences, scaled the
results to $[0,1]$, and applied a threshold of $0.7$ for binarization. Figure
6 displays positions with the highest attention, comparing inferred positions
with mutated positions in the humanized sequence of Certolizumab (see Section
C in the appendix for more results). The model focuses on positions initially
mutated in the parental sequence, recommending fewer mutations than those
found in the humanized sequences. This aligns with observations in (Marks et
al., 2021), suggesting the potential for fewer mutations in humanization.

Figure 4: Comparing SelfPAD-inferred (bottom) to experimental (top) positions
in humanizing Certolizumab’s heavy chain.

Ablation study We performed an ablation study on the 217 immunogenicity
dataset to analyze design choices in FT-SelfPAD. Table 6 shows that
incorporating shift and swap noise, along with re-initializing four pre-
trained model layers before fine-tuning, improves performance. Notably, the
model is less effective when trained from scratch or fine-tuned without random
initialization of its four layers (see first and second rows in Table 6).

Table 6: Ablation study on 217 immunogenicity data. Swap: Swap noise, Shift:
Shift noise, RI: Four pre-trained model layers randomly initialized before
fine-tuning, FS: Training from scratch without pre-training.

FS | RI | Swap | Shift | ROC AUC | PR AUC
---|---|---|---|---|---
+ | NA | + | + | 80.32 | 92.32
NA | - | + | + | 81.90 | 93.33
NA | - | - | - | 80.74 | 91.89
NA | + | - | - | 81.53 | 92.21
NA | + | + | - | 82.37 | 93.54
NA | + | + | + | 82.86 | 93.79

##  5 Conclusion

In this work, we present a framework, SelfPAD, to learn representations from
the Patented Antibody Database, and demonstrate its usefulness by providing a
robust model for humanness prediction - a practical challenge in therapeutic
antibody development. SelfPAD opens new avenues for advancing the
understanding and optimization of antibody sequences, ultimately facilitating
the development of more effective and safe therapeutic interventions. Despite
its advantages, one of the main drawbacks of our approach is that the PAD is a
noisy dataset and further curation is needed to improve the signal-to-noise
ratio. Moreover, the PAD introduces new biases e.g. it biases towards lab
optimised rather than natural sequences. Thus, we can also improve our results
by incorporating data from other sources such as the OAS since it includes a
large number of natural sequences. Finally, although this work focuses on
humanness prediction, we can fine-tune SelfPAD for other downstream tasks such
as predicting developability properties of antibodies.

## References

  * Abhinandan & Martin (2007) Abhinandan, K. and Martin, A. C.  Analyzing the “degree of humanness” of antibody sequences.  _Journal of molecular biology_ , 369(3):852–862, 2007.
  * Caron et al. (2020) Caron, M., Misra, I., Mairal, J., Goyal, P., Bojanowski, P., and Joulin, A.  Unsupervised learning of visual features by contrasting cluster assignments.  _arXiv preprint arXiv:2006.09882_ , 2020.
  * Chen et al. (2020a) Chen, T., Kornblith, S., Norouzi, M., and Hinton, G.  A simple framework for contrastive learning of visual representations.  In _International conference on machine learning_ , pp. 1597–1607. PMLR, 2020a.
  * Chen et al. (2020b) Chen, X., Fan, H., Girshick, R., and He, K.  Improved baselines with momentum contrastive learning.  _arXiv preprint arXiv:2003.04297_ , 2020b.
  * Clavero-Álvarez et al. (2018) Clavero-Álvarez, A., Di Mambro, T., Perez-Gaviro, S., Magnani, M., and Bruscolini, P.  Humanization of antibodies using a statistical inference approach.  _Scientific reports_ , 8(1):14820, 2018.
  * Collobert & Weston (2008) Collobert, R. and Weston, J.  A unified architecture for natural language processing: Deep neural networks with multitask learning.  In _Proceedings of the 25th international conference on Machine learning_ , pp. 160–167, 2008.
  * Conneau et al. (2019) Conneau, A., Khandelwal, K., Goyal, N., Chaudhary, V., Wenzek, G., Guzmán, F., Grave, E., Ott, M., Zettlemoyer, L., and Stoyanov, V.  Unsupervised cross-lingual representation learning at scale.  _arXiv preprint arXiv:1911.02116_ , 2019.
  * Devlin et al. (2018) Devlin, J., Chang, M.-W., Lee, K., and Toutanova, K.  Bert: Pre-training of deep bidirectional transformers for language understanding.  _arXiv preprint arXiv:1810.04805_ , 2018.
  * Ecker et al. (2015) Ecker, D. M., Susana D., J., and Howard L., L.  The therapeutic monoclonal antibody market.  _MAbs._ , 2015.
  * Falcon & Cho (2020) Falcon, W. and Cho, K.  A framework for contrastive self-supervised learning and designing a new approach.  _arXiv preprint arXiv:2009.00104_ , 2020.
  * Gao et al. (2013) Gao, S. H., Huang, K., Tu, H., and Adler, A. S.  Monoclonal antibody humanness score and its applications.  _BMC biotechnology_ , 13(1):1–12, 2013.
  * Grill et al. (2020) Grill, J.-B., Strub, F., Altché, F., Tallec, C., Richemond, P. H., Buchatskaya, E., Doersch, C., Pires, B. A., Guo, Z. D., Azar, M. G., et al.  Bootstrap your own latent: A new approach to self-supervised learning.  _arXiv preprint arXiv:2006.07733_ , 2020.
  * He et al. (2020) He, K., Fan, H., Wu, Y., Xie, S., and Girshick, R.  Momentum contrast for unsupervised visual representation learning.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 9729–9738, 2020.
  * Joulin et al. (2016) Joulin, A., Grave, E., Bojanowski, P., and Mikolov, T.  Bag of tricks for efficient text classification.  _arXiv preprint arXiv:1607.01759_ , 2016.
  * Krawczyk et al. (2021) Krawczyk, K., Buchanan, A., and Marcatili, P.  Data mining patented antibody sequences.  In _MAbs_ , volume 13, pp. 1892366. Taylor & Francis, 2021.
  * Liu et al. (2019) Liu, Y., Ott, M., Goyal, N., Du, J., Joshi, M., Chen, D., Levy, O., Lewis, M., Zettlemoyer, L., and Stoyanov, V.  Roberta: A robustly optimized bert pretraining approach.  _arXiv preprint arXiv:1907.11692_ , 2019.
  * Loshchilov & Hutter (2017) Loshchilov, I. and Hutter, F.  Decoupled weight decay regularization.  _arXiv preprint arXiv:1711.05101_ , 2017.
  * Marks et al. (2021) Marks, C., Hummer, A. M., Chin, M., and Deane, C. M.  Humanization of antibodies using a machine learning approach on large-scale repertoire data.  _Bioinformatics_ , 37(22):4041–4047, 2021.
  * Mikolov et al. (2013) Mikolov, T., Chen, K., Corrado, G., and Dean, J.  Efficient estimation of word representations in vector space.  _arXiv preprint arXiv:1301.3781_ , 2013.
  * Olsen et al. (2022a) Olsen, T. H., Boyles, F., and Deane, C. M.  Observed antibody space: A diverse database of cleaned, annotated, and translated unpaired and paired antibody sequences.  _Protein Science_ , 31(1):141–146, 2022a.
  * Olsen et al. (2022b) Olsen, T. H., Moal, I. H., and Deane, C. M.  Ablang: an antibody language model for completing antibody sequences.  _Bioinformatics Advances_ , 2(1):vbac046, 2022b.
  * Oord et al. (2018) Oord, A. v. d., Li, Y., and Vinyals, O.  Representation learning with contrastive predictive coding.  _arXiv preprint arXiv:1807.03748_ , 2018.
  * Paszke et al. (2019) Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., and Chintala, S.  Pytorch: An imperative style, high-performance deep learning library.  In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems 32_ , pp. 8024–8035. Curran Associates, Inc., 2019.  URL http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf.
  * Pennington et al. (2014) Pennington, J., Socher, R., and Manning, C. D.  Glove: Global vectors for word representation.  In _Proceedings of the 2014 conference on empirical methods in natural language processing (EMNLP)_ , pp. 1532–1543, 2014.
  * Poiron et al. (2010) Poiron, C., Wu, Y., Ginestoux, C., Ehrenmann, F., Duroux, P., and Lefranc, M.  Imgt/mab-db: the imgt® database for therapeutic monoclonal antibodies.  _Poster no101_ , 11:382, 2010.
  * Prihoda et al. (2022) Prihoda, D., Maamary, J., Waight, A., Juan, V., Fayadat-Dilman, L., Svozil, D., and Bitton, D. A.  Biophi: A platform for antibody design, humanization, and humanness evaluation based on natural antibody repertoires and deep learning.  In _MAbs_ , volume 14, pp. 2020203. Taylor & Francis, 2022.
  * Radford et al. (2019) Radford, A., Wu, J., Child, R., Luan, D., Amodei, D., Sutskever, I., et al.  Language models are unsupervised multitask learners.  _OpenAI blog_ , 1(8):9, 2019.
  * Ramon et al. (2024) Ramon, A., Ali, M., Atkinson, M., Saturnino, A., Didi, K., Visentin, C., Ricagno, S., Xu, X., Greenig, M., and Sormanni, P.  Assessing antibody and nanobody nativeness for hit selection and humanization with abnativ.  _Nature Machine Intelligence_ , Jan 2024.  ISSN 2522-5839.  doi: 10.1038/s42256-023-00778-3.  URL https://doi.org/10.1038/s42256-023-00778-3.
  * Ruffolo et al. (2021) Ruffolo, J. A., Gray, J. J., and Sulam, J.  Deciphering antibody affinity maturation with language models and weakly supervised learning.  _arXiv preprint arXiv:2112.07782_ , 2021.
  * Shuai et al. (2023) Shuai, R. W., Ruffolo, J. A., and Gray, J. J.  Iglm: Infilling language modeling for antibody sequence design.  _Cell Systems_ , 14(11):979–989, 2023.
  * Ucar et al. (2021) Ucar, T., Hajiramezanali, E., and Edwards, L.  Subtab: Subsetting features of tabular data for self-supervised representation learning.  _Advances in Neural Information Processing Systems_ , 34:18853–18865, 2021.
  * Van Den Oord et al. (2017) Van Den Oord, A., Vinyals, O., et al.  Neural discrete representation learning.  _Advances in neural information processing systems_ , 30, 2017.
  * Vandivort et al. (2020) Vandivort, T. C., Horton, D. B., and Johnson, S. B.  Regulatory and strategic considerations for addressing immunogenicity and related responses in biopharmaceutical development programs.  _Journal of Clinical and Translational Science_ , 4(6):547–555, 2020.

##  Appendix A Amino Acid Features

Table 7: 18 biophysical features of amino acids used in SelfPAD.

Hydropathy Chemical Properties Volume Charge Hydrogen Polarity Hydrophobic
Neutral Hyrophilic Aliphatic Sulfur Hydroxyl Amide Basic Acidic Aromatic
Positive Uncharged Negative Donor Acceptor Polar Nonpolar I 4.5 0 0 1 0 0 0 0
0 0 166.7 0 1 0 0 0 0 1 V 4.2 0 0 1 0 0 0 0 0 1 140 0 1 0 0 0 0 1 L 3.8 0 0 1
0 0 0 0 0 0 166.7 0 1 0 0 0 0 1 F 2.8 0 0 0 0 0 0 0 0 1 189.9 0 1 0 0 0 0 1 C
2.5 0 0 0 1 0 0 0 0 0 108.5 0 1 0 0 0 0 1 M 1.9 0 0 0 1 0 0 0 0 0 162.9 0 1 0
0 0 0 1 A 1.8 0 0 1 0 0 0 0 0 0 88.6 0 1 0 0 0 0 1 W* -0.9 0 0 0 0 0 0 0 0 1
227.8 0 1 0 1 0 0 1 G 0 -0.4 0 0 0 0 0 0 0 0 60.1 0 1 0 0 0 0 1 T 0 -0.7 0 0 0
1 0 0 0 0 116.1 0 1 0 1 1 1 0 S 0 -0.8 0 0 0 1 0 0 0 0 89 0 1 0 1 1 1 0 Y 0
-1.3 0 0 0 0 0 0 0 0 193.6 0 1 0 1 1 1 0 P 0 -1.6 0 0 0 0 0 0 0 0 112.7 0 1 0
0 0 0 1 H 0 -3.2 0 0 0 0 0 1 0 0 153.2 1 0 0 1 1 1 0 N 0 0 -3.5 0 0 0 1 0 0 0
114.1 0 1 0 1 1 1 0 D 0 0 -3.5 0 0 0 0 0 1 0 111.1 0 0 1 0 1 1 0 Q 0 0 -3.5 0
0 0 1 0 0 0 143.8 0 1 0 1 1 1 0 E 0 0 -3.5 0 0 0 0 0 1 0 138.4 0 0 1 0 1 1 0 K
0 0 -3.9 0 0 0 0 1 0 0 168.6 1 0 0 1 0 1 0 R 0 0 -4.5 0 0 0 0 1 0 0 173.4 1 0
0 1 0 1 0

Figure 5: Training loss during pre-training of SelfPAD.

##  Appendix B Implementation and resources

We implemented our work using PyTorch (Paszke et al., 2019). AdamW optimizer
(Loshchilov & Hutter, 2017) with $betas=(0.9,0.999)$ and $eps=1e-07$ is used
for all of our experiments. We used a compute cluster consisting of A10G GPUs
throughout this work. For pre-training, we did a hyper-parameter sweep over
the parameters such as number of layers of the encoder, number of attention
heads, embedding dimension, temperature constant for contrastive loss,
learning rate, total epochs etc. and chose the parameters based on the best
convergence on the training set. For fine-tuning phase, we did a hyper-
parameter search using a validation set.

Table 8 lists hyperparameters used for pre-training and fine-tuning of the
model during this work while Table 9 shows the architecture used for the
encoder. MLP consists of one block with a residual connection, where the block
has a linear layer (128x512) followed by GELU activation, batch normalisation,
dropout (p=0.3) and another linear layer (512x128). A final linear layer
(128x2) with the softmax activation is used to make a binary classification.

Table 8: Hyperparameters: Hyperparameters used for both pre-training and fine-
tuning stages. Abbreviations are; LR: Learning Rate, P: Dropout rate, LS:
Label Smoothing used for cross-entropy loss, specifying the amount of
smoothing when computing the loss, BS: Batch Size, temp.: Temperature
parameter used for contrastive loss, $\bm{P}_{swap}:$ Percentage of embedding
dimension used for swap-noise, $\bm{R}_{shift}:$ Range used to sample a number
that is used to shift the position of sequence to add shift-noise.

| epoch | bs | lr | p | ls | temp. | $\bm{P}_{swap}$ | $\bm{N}_{shift}$
---|---|---|---|---|---|---|---|---
Pre-training | 1000 | 100 | 1e-3 | 0.2 | NA | 0.1 | NA | NA
Fine-tuning | 25 | 512 | [1e-4, 1e-3] | 0.3 | 0.5 | NA | 0.2 | $[0,10]$

Table 9: Encoder architecture: Abbreviations are; $\bm{N}_{layers}:$ Number of
layers, $\bm{N}_{heads}:$ Number of attention heads in each layer,
$\bm{d}_{emb}:$ Dimension of embedding in each layer, $\bm{d}_{head}:$ Feature
dimension of each attention head.

| $\bm{N}_{layers}$ | $\bm{N}_{heads}$ | $\bm{d}_{emb}$ | $\bm{d}_{head}$
---|---|---|---|---
Encoder | 4 | 4 | 32 | 16

##  Appendix C Interpretability

Figure 6: Comparing SelfPAD-inferred positions to experimental data where the
changes are made to humanize the heavy chain of the parental sequences of
Herceptin (top), Bevacizumab (middle) and Pinatuzumab (bottom).


# A closed-form multigrid smoothing factor for an additive Vanka-type smoother
applied to the Poisson equation

Chen Greif Department of Computer Science, The University of British Columbia,
Vancouver, BC, V6T 1Z4, Canada. The work of first author was supported in part
by a Discovery Grant of the Natural Sciences and Engineering Research Council
of Canada<EMAIL_ADDRESS><EMAIL_ADDRESS> Yunhui He11footnotemark: 1

###### Abstract

We consider an additive Vanka-type smoother for the Poisson equation
discretized by the standard finite difference centered scheme. Using local
Fourier analysis, we derive analytical formulas for the optimal smoothing
factors for two types of smoothers, called vertex-wise and element-wise Vanka
smoothers, and present the corresponding stencils. Interestingly, in one
dimension the element-wise Vanka smoother is equivalent to the scaled mass
operator obtained from the linear finite element method, and in two dimensions
the element-wise Vanka smoother is equivalent to the scaled mass operator
discretized by bilinear finite element method plus a scaled identity operator.
Based on these discoveries, the mass matrix obtained from finite element
method can be used as an approximation to the inverse of the Laplacian, and
the resulting mass-based relaxation scheme features small smoothing factors in
one, two, and three dimensions. Advantages of the mass operator are that the
operator is sparse and well conditioned, and the computational cost of the
relaxation scheme is only one matrix-vector product; there is no need to
compute the inverse of a matrix. These findings may help better understand the
efficiency of additive Vanka smoothers and develop fast solvers for numerical
solutions of partial differential equations.

Keywords. local Fourier analysis, multigrid, smoothing factor, additive Vanka-
type smoother, mass matrix, finite difference method

##  1 Introduction

Consider the Poisson equation in one dimension (1D) or two dimensions (2D):

$-\Delta u=f,$ (1)

where $u=u(x)$ and $f=f(x)$ in 1D or $u=u(x,y)$ and $f=f(x,y)$ in 2D. The
function $f$ is assumed to be sufficiently smooth so that finite-difference
discretizations of $u$ provide effective approximations of the solution. The
differential operator $\Delta$ stands for the Laplacian:
$\Delta\equiv\frac{d^{2}}{dx^{2}}$ in 1D and
$\Delta\equiv\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial
y^{2}}$ in 2D. We assume a uniform mesh discretization with meshsize $h$.

We apply the standard three- and five-point finite difference discretizations
for the Laplacian; the corresponding stencils are given by

$A_{h}=\frac{1}{h^{2}}\begin{bmatrix}-1&2&-1\end{bmatrix}$ (2)

and

$A_{h}=\frac{1}{h^{2}}\begin{bmatrix}&-1&\\\ -1&4&-1\\\ &-1&\end{bmatrix},$
(3)

respectively.

Let us denote the corresponding linear system by

$A_{h}u_{h}=b_{h}.$ (4)

The numerical solution of (4) is one of the most extensively explored topics
in the numerical linear algebra literature. The discrete Laplacian $A_{h}$ is
a symmetric positive definite M-matrix, its eigenvalues are explicitly known,
and it is used as a primary benchmark problem for the development of fast
solvers. When the problem is large, iterative solvers that allow a high level
of parallelism are often preferred.

One of the most efficient methods for solving (4) is multigrid [14, 18]. The
choice of multigrid components, such as a relaxation scheme as a smoother,
plays an important role in the design of fast methods.

The choice of additive Vanka as relaxation scheme, which is suitable for
parallel computing, has recently drawn considerable attention. Vanka-type
smoothers have been applied to the Navier–Stokes equations [10, 12, 16], the
Poisson equation using continuous and discontinuous finite elements methods
[9], the Stokes equations with finite element methods [6, 17], poroelasticity
equations [2] in monolithic multigrid, and other problems. A restricted
additive Vanka for Stokes using $Q_{1}-Q_{1}$ discretizations is presented in
[13], which shows its competitiveness with the multiplicative Vanka smoother.
Nonoverlapping block smoothing using different patches has been discussed in
[3] for the Stokes equations discretized by the marker-and-cell scheme.
Muliplicative Vanka smoothers in combination with multigrid methods are
discussed in [4, 5, 7]. Vanka-type relaxation has been used in many contexts,
and for more details, we refer the reader to [1, 10, 11].

In the literature, Vanka-type relaxation schemes demonstrate their high
efficiency in a multigrid setting, but there seems no theoretical analysis for
the convergence speed even for the simple Poisson equation. In this work we
take steps towards closing this gap by considering the additive Vanka
relaxation for the Poisson equation, and exploring stencils for the Vanka
patches. We derive analytical optimal smoothing factors for two types of
additive Vanka patches used in [9] for hybridized and embedded discontinuous
Galerkin methods. Moreover, we also find the corresponding stencils for the
Vanka operator, and show that they are closely related to the scaled mass
matrix obtained from the finite element method. Based on this discovery, we
propose the mass-based relaxation scheme, which yields rapid convergence. This
mass-based relaxation is very simple: the computation cost is only matrix-
vector product and there is no need to solve the subproblems needed in an
additive Vanka setting. Another advantage of the mass matrices obtained from
(bi)linear elements is sparsity.

Solvers for the Poisson equation often form the first step for designing fast
solvers for more complex problems, such as the Stokes equations, and Navier-
Stokes equations. We therefore believe that the findings in this work can give
some hints in future for designing fast numerical methods for these complex
problems.

The remainder of this paper is organized as follows. In Section 2 we introduce
the two types of additive Vanka smoothers for the Poisson equation. In Section
3, we present our theoretical analysis of optimal smoothing factors in one and
two dimensions. Based on our analysis we also propose a mass-based smoother
for the three-dimensional problem, where the mass matrix is obtained from the
trilinear finite element method. In Section 4, we numerically validate our
analytical observations and present an LFA two-grid convergence factor.
Finally, in Section 5 we discuss our findings and draw some conclusions.

##  2 Vanka-type smoother

We are interested in exploring the structure of an additive Vanka-type
smoother for solving the linear system (4) using multigrid. In general, this
type can be thought of as related to the family of block Jacobi smoothers,
which are suitable for parallel computation and are typically highly efficient
within the context of multigrid smoothing.

Let the degrees of freedom (DoFs) of $u_{h}$ be the set $\Upsilon$ such that
$\Upsilon=\bigcup_{i=1}^{N}\Upsilon_{i}$. $V_{i}$ is a restriction operator
that maps the vector $u_{h}$ onto the vector in $\Upsilon_{i}$. Define

$A_{i}=V_{i}A_{h}V^{T}_{i}.$

Then, we update current approximation $u^{j}$ by a single Vanka relaxation
given by:
For $i=1,\cdots,N$,

$A_{i}\delta u_{i}=V_{i}(b_{h}-A_{h}u^{j})$

and

$u^{j+1}=u^{j}+\sum_{i=1}^{N}V_{i}^{T}W_{i}\delta u_{i}.$

A single iteration of Vanka can be represented as

$M=\sum_{i=1}^{N}V_{i}^{T}W_{i}A_{i}^{-1}V_{i},$ (5)

where the weighting matrix $W=W_{i}$ is given by the natural weights of the
overlapping block decomposition. Each diagonal entry is equal to the
reciprocal of the number of patches that the corresponding degree of freedom
appears in. We refer to $M$ the Vanka operator.

For a single additive Vanka relaxation process, the relaxation error operator
is given by

$S=I-\omega MA_{h}.$ (6)

A key factor is the choice of the patch, that is, the $\\{\Upsilon_{i}\\}$.
Here, following [9], we consider two patches, shown in Figure 1. We refer to
the left patch in Figure 1 as an element-wise patch and the right one as a
vertex-wise patch. We denote the corresponding relaxation error operators in
(6) as $S_{e}$ and $S_{v}$, respectively. The collection of circles indicate
the number of DoFs in one patch $\Upsilon_{i}$. This means that the resulting
subproblem is associated with a small matrix $A_{i}$ whose size is $4\times 4$
or $5\times 5$. In the remainder of this work, for simplicity and clarity we
use subscripts $e$ and $v$ to distinguish between the corresponding operators
for element-wise Vanka and vertex-wise Vanka, respectively.

$\bigcirc$$\bigcirc$$\bigcirc$$\bigcirc$

$\bigcirc$$\bigcirc$$\bigcirc$$\bigcirc$$\bigcirc$

Figure 1: Left: element-wise patch. Right: vertex-wise patch.

##  3 Local Fourier analysis

Local Fourier analysis (LFA) [15, 19] is a useful tool for predicting and
analyzing the convergence behavior of multigrid and other numerical
algorithms. In this section, we employ LFA to examine the spectrum or spectral
radius of the underlying operators to better understand the proposed Vanka-
type smoother.

When applying LFA to multigrid methods, high and low frequencies [15] for
standard coarsening ($H=2h$) are introduced:

$\boldsymbol{\theta}\in T^{{\rm
low}}=\left[-\frac{\pi}{2},\frac{\pi}{2}\right)^{d},\quad\boldsymbol{\theta}\in
T^{{\rm
high}}=\displaystyle\left[-\frac{\pi}{2},\frac{3\pi}{2}\right)^{d}\bigg{\backslash}\left[-\frac{\pi}{2},\frac{\pi}{2}\right)^{d},$

where $d$ is the dimension of the underlying problem.

To quantitatively assess the performance of multigrid, LFA-predicted two-grid
convergence and smoothing factors are used. In many cases the LFA smoothing
factor is used due to its simplicity. It works under the assumption that the
smoothing process reduces the high frequencies and leaves the low frequencies
unchanged. The LFA smoothing factor typically offers a rather sharp bound on
the actual two-grid performance.

######  Definition 3.1.

Let $S$ be the relaxation error operator. Then, the corresponding LFA
smoothing factor for $\mathcal{S}$ is given by

$\mu(\omega)=\max_{\boldsymbol{\theta}\in T^{{\rm
high}}}\Big{\\{}\rho(\widetilde{S}(\boldsymbol{\theta},\omega))\,\,\Big{\\}},$
(7)

where $\widetilde{S}(\boldsymbol{\theta},\omega)$ is the symbol of $S$,
$\omega$ is algorithmic parameter, and
$\rho(\widetilde{S}(\boldsymbol{\theta},\omega))$ denotes the spectral radius
of the matrix $\widetilde{S}(\boldsymbol{\theta},\omega)$.

Note that the LFA smoothing factor $\mu$ is a function of $\omega$. Often, one
can minimize (7) with respect to $\omega$ to obtain fast convergence speed. We
define the optimal smoothing factor as

$\mu_{\rm opt}=\min_{\omega}\mu.$ (8)

In this work, the symbol $\widetilde{S}$ for the Laplacian considered is a
scalar, so the spectral radius is reduced to the maximum of a scalar function.
In the following, we use LFA to identify the optimal smoothing factor for the
additive Vanka-type relaxation schemes and explore the structure of the Vanka
operator $M$ defined in (5). Before providing our detailed analysis of the
smoothing factor for different relaxation schemes, we summarize our results in
Table 1. The table provides a review of quantitative results of additive Vanka
smoothers, and for comparison we include results for the standard point-wise
damped Jacobi smoother.

Table 1: The optimal LFA smoothing factor in 1D and 2D Smoother | Jacobi | AS-e | AS-v
---|---|---|---
1D
$\omega_{\rm opt}$ | 2/3 | 12/17 | 81/104
$\mu_{\rm opt}$ | 0.333 | 0.059 | 0.039
2D
$\omega_{\rm opt}$ | 4/5 | 24/25 | 20/23
$\mu_{\rm opt}$ | 0.600 | 0.280 | 0.391

Note that a general form of the symbol of additive Vanka operator for the
Stokes equations has been discussed in [6], which gives
$\widetilde{M}=\widetilde{V}^{T}\widetilde{W}\Phi^{T}A_{i}^{-1}\Phi\widetilde{V}$,
where $\Phi$ is called the relative Fourier matrix. Here, we can directly
apply the formula of $\widetilde{M}$ to our additive Vanka operator; see [6].

In the analysis that follows, we will use $\iota$ to denote the imaginary
scalar satisfying

$\iota^{2}=-1.$

###  3.1 Symbols of Vanka patches in 1D

In this subsection, we first consider the analytical symbol of the element-
wise patch, then the vertex-wise patch for the Laplacian in 1D. We discuss the
optimal smoothing factor for each case and derive the corresponding stencil
for the Vanka operator.

####  3.1.1 Element-wise Vanka patch in 1D

It can easily be shown that the symbol of $A_{h}$, see (2), is given by

$\widetilde{A}_{h}=\frac{1}{h^{2}}2(1-\cos\theta).$ (9)

Moreover, for the element-wise patch the subproblem matrix is

$A_{i}=\frac{1}{h^{2}}\begin{pmatrix}2&-1\\\ -1&2\end{pmatrix}.$

Following [6], the relative Fourier matrix $\Phi$ is

$\Phi=\begin{pmatrix}1&0\\\ 0&e^{\iota\theta}\end{pmatrix}.$

Then, the symbol of $M_{e}$ is given by
$\widetilde{M}_{e}=\widetilde{V}^{T}\widetilde{W}\Phi^{T}A_{i}^{-1}\Phi\widetilde{V}$,
where

$\widetilde{V}=\begin{pmatrix}1\\\
1\end{pmatrix},\quad\widetilde{W}=\frac{1}{2}\begin{pmatrix}1&0\\\
0&1\end{pmatrix},\quad A_{i}^{-1}=\frac{h^{2}}{3}\begin{pmatrix}2&1\\\
1&2\end{pmatrix}.$

Based on the above formulas, we obtain

$\widetilde{M}_{e}=\frac{h^{2}}{6}(4+e^{\iota\theta}+e^{-\iota\theta}).$ (10)

Formula (10) indicates that the element-wise Vanka patch corresponds to the
stencil

$M_{e}=\frac{h^{2}}{6}\begin{bmatrix}1&4&1\end{bmatrix}.$ (11)

Recall that the mass stencil in 1D using linear finite elements is given by

$M_{fe}=\frac{h}{6}\begin{bmatrix}1&4&1\end{bmatrix}.$ (12)

This means that the element-wise Vanka operator is equivalent to a scaled mass
matrix obtained from the linear finite element method, $M_{e}=hM_{fe}$.

Next, we give the optimal smoothing factor for the element-wise Vanka
relaxation scheme.

######  Theorem 3.1.

The optimal smoothing factor of $S_{e}$ for the vertex-wise Vanka in 1D is

$\mu_{e,\rm opt}=\min_{\omega}\max_{\theta\in T^{\rm
high}}{|1-\omega\widetilde{M}_{e}\widetilde{A}_{h}|}=\frac{1}{17}\approx
0.059,$ (13)

where the minimum is uniquely achieved at $\omega=\omega_{\rm
opt}=\frac{12}{17}\approx 0.706$.

###### Proof.

When $\theta\in T^{\rm high}=[\frac{\pi}{2},\frac{3\pi}{2}]$,

$\widetilde{M}_{e}\widetilde{A}_{h}=\frac{2}{3}(2-\cos\theta-\cos^{2}\theta)\in\left[\frac{4}{3},\frac{3}{2}\right].$

Thus,

$\mu(\omega)=\max{|1-\omega\widetilde{M}_{e}\widetilde{A}_{h}|}=\max\left\\{\left|1-\omega\frac{4}{3}\right|,\left|1-\omega\frac{3}{2}\right|\right\\}.$

To minimize $\mu(\omega)$, we require

$\left|1-\omega\frac{4}{3}\right|=\left|1-\omega\frac{3}{2}\right|,$

which gives $\omega=\frac{2}{4/3+3/2}=\frac{12}{17}$. Then, $\mu_{\rm
opt}=1-\frac{2}{4/3+3/2}\frac{4}{3}=\frac{1}{17}$. ∎

It is well known that the optimal smoothing factor for damped Jacobi
relaxation for the Laplacian in 1D (with $\omega=\frac{2}{3}$) is
$\frac{1}{3}\approx 0.333\gg 0.059$. This suggests that using the additive
Vanka smoother for multigrid achieves much faster convergence.

####  3.1.2 Vertex-wise Vanka patches in 1D

We now consider the vertex-wise patch. The subproblem matrix is

$A_{i}=\frac{1}{h^{2}}\begin{pmatrix}2&-1&0\\\ -1&2&-1\\\
0&-1&2\end{pmatrix}.$

Again, following [6], the relative Fourier matrix $\Phi$ is

$\Phi=\begin{pmatrix}e^{-\iota\theta}&0&0\\\ 0&1&0\\\
0&0&e^{\iota\theta}\end{pmatrix}.$

Then, the symbol of $M_{v}$ is given by
$\widetilde{M}_{v}=\widetilde{V}^{T}\widetilde{W}\Phi^{T}A_{i}^{-1}\Phi\widetilde{V}$,
where

$\widetilde{V}=\begin{pmatrix}1\\\ 1\\\
1\end{pmatrix},\quad\widetilde{W}=\frac{1}{3}\begin{pmatrix}1&0&0\\\ 0&1&0\\\
0&0&1\end{pmatrix},\quad A_{i}^{-1}=\frac{h^{2}}{4}\begin{pmatrix}3&2&1\\\
2&4&2\\\ 1&2&3\end{pmatrix}.$

Using the above formulas, we have

$\widetilde{M}=\widetilde{V}^{T}\widetilde{W}\Phi^{T}A_{i}^{-1}\Phi\widetilde{V}=\frac{h^{2}}{12}(10+4e^{\iota\theta}+4e^{-\iota\theta}+e^{\iota
2\theta}+4e^{-\iota 2\theta}).$ (14)

Based on (14), the stencil of $M_{v}$ is

$M_{v}=\frac{h^{2}}{12}\begin{bmatrix}1&4&10&4&1\end{bmatrix}.$ (15)

Compared with (11), the vertex-wise Vanka uses a wider stencil.

######  Theorem 3.2.

The optimal smoothing factor of $S_{v}$ for the vertex-wise Vanka in 1D is

$\mu_{v,\rm opt}=\min_{\omega}\max_{\theta\in T^{\rm
high}}{|1-\omega\widetilde{M}_{v}\widetilde{A}_{h}|}=\frac{1}{26}\approx
0.039,$ (16)

where the minimum is uniquely achieved at $\omega=\omega_{\rm
opt}=\frac{81}{104}\approx 0.779$.

###### Proof.

From (9) and (14), we have

$\displaystyle\widetilde{M}_{v}\widetilde{A}_{h}$
$\displaystyle=\frac{1}{3}(1-\cos\theta)\left(5+4\cos\theta+\cos(2\theta)\right),$
$\displaystyle=\frac{2}{3}(1-\cos\theta)(2+2\cos\theta+\cos^{2}\theta),$
$\displaystyle=\frac{2}{3}(1-\cos\theta)\left((\cos\theta+1)^{2}+1\right).$

Note that when $\theta\in T^{\rm
high}=[\frac{\pi}{2},\frac{3\pi}{2}],\cos\theta\in[-1,0]$. Let
$g(x)=\frac{2}{3}(1-x)\left((x+1)^{2}+1\right)$, where
$x=\cos\theta\in[-1,0]$. To identify the range of $g(x)$, we first compute its
derivative:

$g^{\prime}(x)=\frac{-2}{3}x(3x+2)<0,\quad\text{for}\,\,x\in[-1,0].$

It follows that

$\left(\frac{10}{9}\right)^{2}=g(-2/3)\leq g(x)\leq g(0)=g(1)=\frac{4}{3}.$

That is,
$\widetilde{M}_{v}\widetilde{A}_{h}\in[(\frac{10}{9})^{2},\frac{4}{3}]$ for
$\theta\in T^{\rm high}$ . To minimize $\mu(\omega)$, we require
$\omega=\frac{2}{4/3+(10/9)^{2}}=\frac{81}{104}$. Then, it follows that
$\mu_{\rm opt}=1-\frac{2}{4/3+(10/9)^{2}}(\frac{10}{9})^{2}=\frac{1}{26}.$

∎

Again, the optimal smoothing factor for vertex-wise Vanka is significantly
smaller (and hence better) than that of the damped Jacobi relaxation scheme.

###  3.2 Symbols of Vanka patches in 2D

Similarly to previous subsection, we first consider the analytical symbol for
the element-wise patch, then for the vertex-wise patch for the Laplacian in
2D.

####  3.2.1 Element-wise Vanka patch in 2D

The symbol of Laplace operator discretized by five-point stencil, see (3), is

$\widetilde{A}_{h}=\frac{1}{h^{2}}(4-e^{\iota\theta_{1}}-e^{\iota\theta_{2}}-e^{-\iota\theta_{1}}-e^{-\iota\theta_{2}})=\frac{1}{h^{2}}\big{(}4-2\cos\theta_{1}-2\cos\theta_{2}\big{)}.$
(17)

For the vertex-wise patch, following [6], the relative Fourier matrix $\Phi$
is

$\Phi=\begin{pmatrix}1&0&0&0\\\ 0&e^{\iota\theta_{1}}&0&0\\\
0&0&e^{\iota\theta_{2}}&0\\\
0&0&0&e^{\iota(\theta_{1}+\theta_{2})}\end{pmatrix}.$

Then, the symbol of $M_{e}$ is given by
$\widetilde{M}_{e}=\widetilde{V}^{T}\widetilde{W}\Phi^{T}A_{i}^{-1}\Phi\widetilde{V}$,
where

$\widetilde{W}=\frac{1}{4}\begin{pmatrix}1&0&0&0\\\ 0&1&0&0\\\ 0&0&1&0\\\
0&0&0&1\end{pmatrix},\quad\widetilde{V}=\begin{pmatrix}1\\\ 1\\\ 1\\\
1\end{pmatrix},\quad A_{i}=\frac{1}{h^{2}}\begin{pmatrix}4&-1&-1&0\\\
-1&4&0&-1\\\ -1&0&4&-1\\\ 0&-1&-1&4\end{pmatrix}.$

We have

$A_{i}^{-1}=h^{2}\begin{pmatrix}7/24&1/12&1/12&1/24\\\ 1/12&7/24&1/24&1/12\\\
1/12&1/24&7/24&1/12\\\ 1/24&1/12&1/12&7/24\end{pmatrix},$

and it follows that

$\displaystyle\widetilde{M}_{e}$
$\displaystyle=\widetilde{V}^{T}\widetilde{W}\Phi^{T}A_{i}^{-1}\Phi\widetilde{V},$
$\displaystyle=\frac{h^{2}}{96}\big{(}28+4(e^{\iota\theta_{1}}+e^{-\iota\theta_{1}}+e^{\iota\theta_{2}}+e^{-\iota\theta_{2}})+e^{\iota(\theta_{1}+\theta_{2})}+e^{-\iota(\theta_{1}+\theta_{2})}+(e^{\iota\theta_{1}}e^{-\iota\theta_{2}}+e^{-\iota\theta_{1}}e^{\iota\theta_{2}})\big{)},$
$\displaystyle=\frac{h^{2}}{24}\big{(}7+2(\cos\theta_{1}+\cos\theta_{2})+\cos\theta_{1}\cos\theta_{2}\big{)}.$
(18)

Next, we give the optimal smoothing factor for the element-wise Vanka
relaxation in 2D.

######  Theorem 3.3.

The optimal smoothing factor of $S_{e}$ for the element-wise Vanka relaxation
in 2D is given by

$\mu_{e,\rm opt}=\min_{\omega}\max_{\boldsymbol{\theta}\in T^{{\rm
high}}}{|1-\omega\widetilde{M}_{e}\widetilde{A}_{h}|}=\frac{7}{25}\approx
0.280,$ (19)

where the minimum is uniquely achieved at $\omega=\omega_{\rm
opt}=\frac{24}{25}\approx 0.960$.

###### Proof.

From (17) and (18), we have

$\displaystyle\widetilde{M}_{e}\widetilde{A}_{h}$
$\displaystyle=\frac{1}{24}\big{(}7+2(\cos\theta_{1}+\cos\theta_{2})+\cos\theta_{1}\cos\theta_{2}\big{)}\big{(}4-2\cos\theta_{1}-2\cos\theta_{2}\big{)},$
$\displaystyle=\frac{1}{12}\big{(}7+2(\cos\theta_{1}+\cos\theta_{2})+\cos\theta_{1}\cos\theta_{2}\big{)}\big{(}2-(\cos\theta_{1}+\cos\theta_{2})\big{)},$
$\displaystyle=\frac{1}{12}\big{(}7+2(\eta_{1}+\eta_{2})+\eta_{1}\eta_{2}\big{)}\big{(}2-\eta_{1}-\eta_{2}\big{)},$

where $\eta_{1}=\cos\theta_{1},\eta_{2}=\cos\theta_{2}$. Let

$g(\eta_{1},\eta_{2})=\frac{1}{12}\big{(}7+2(\eta_{1}+\eta_{2})+\eta_{1}\eta_{2}\big{)}\big{(}2-\eta_{1}-\eta_{2}\big{)}.$

We have $(\eta_{1},\eta_{2})\in[-1,-1]^{2}=:\Omega$ and $g(\eta_{1},\eta_{2})$
is continuous in $\Omega$. By the Extreme Value Theorem, $g$ achieves its
extremal values at the boundary of $\Omega$ or its derivatives are zeros. We
first consider the derivatives,

$\displaystyle g_{\eta_{1}}$
$\displaystyle=\frac{1}{12}(-3-4\eta_{1}-2\eta_{1}\eta_{2}-\eta_{2}^{2}-2\eta_{2}),$
$\displaystyle g_{\eta_{2}}$
$\displaystyle=\frac{1}{12}(-3-4\eta_{2}-2\eta_{1}\eta_{2}-\eta_{2}^{2}-2\eta_{1}).$

Solving $g_{\eta_{1}}=g_{\eta_{1}}=0$ gives $\eta_{1}=\eta_{2}=-1$. Thus,
$g(-1,-1)=\frac{4}{3}$ is a possible global extreme value.

Next, we compute the extremal values of $g(\eta_{1},\eta_{2})$ at the boundary
$\Omega$. Due to the symmetric of $g$, we only need to consider the following
two cases.

  * •

Case 1: $\eta_{1}=-1$ and $\eta_{2}\in[-1,1]$. We have

$g(-1,\eta_{2})=\frac{1}{12}(5+\eta_{2})(3-\eta_{2}).$

Thus,

$1=g(-1,1)\leq g(-1,\eta_{2})\leq g(-1,-1)=\frac{4}{3}.$

  * •

Case 2: $\eta_{1}=1$ and $\eta_{2}\in[-1,1]$. We have

$g(1,\eta_{2})=\frac{1}{4}(3+\eta_{2})(1-\eta_{2}).$

Thus,

$0=g(1,1)\leq g(1,\eta_{2})\leq g(1,-1)=1.$

If we restrict $(\theta_{1},\theta_{2})\in T^{\rm high}$, we find that the
maximum and minimum of
$\widetilde{M}_{e}\widetilde{A}_{h}=g(\eta_{1},\eta_{2})$ are given by

$g(-1,-1)=\frac{4}{3},\quad g(1,0)=\frac{3}{4},$ (20)

respectively. It follows that $\omega_{\rm
opt}=\frac{2}{3/4+4/3}=\frac{24}{25}$ and $\mu_{e,\rm
opt}=1-\frac{2}{4/3+3/4}\frac{3}{4}=\frac{7}{25}$. ∎

It is well known that the optimal smoothing factor for damped Jacobi
relaxation for the Laplacian in 2D is $\frac{3}{5}$ with $\omega=\frac{4}{5}$
[15]. This suggests that using the additive Vanka smoother for multigrid
method, convergence is faster compared to the damped Jacobi relaxation scheme.

Based on the symbol of $M_{e}$, we obtain the stencil of $M_{e}$,

$M_{e}=\frac{h^{2}}{96}\begin{bmatrix}1&4&1\\\ 4&28&4\\\ 1&4&1\end{bmatrix}.$
(21)

Recall that the mass matrix stencil using bilinear finite elements is

$M_{fe}=\frac{h^{2}}{36}\begin{bmatrix}1&4&1\\\ 4&16&4\\\ 1&4&1\end{bmatrix}.$
(22)

Now, we can make a connection between $M_{e}$ and $M_{fe}$:

$M_{e}=\frac{3}{8}M_{fe}+\frac{h^{2}}{8}\mathcal{I},$ (23)

where

$\mathcal{I}=\begin{bmatrix}0&0&0\\\ 0&1&0\\\ 0&0&0\end{bmatrix}.$ (24)

The relationship (23) is interesting, and suggests that the mass matrix
obtained from bilinear elements might be a good approximation to the inverse
of $A_{h}$. Let us, then, move to consider the mass matrix (22) as an
approximation to the inverse of $A_{h}$.

######  Theorem 3.4.

Given the mass stencil $M_{fe}$ in (22) and the relaxation scheme
$S_{fe}=I-\omega M_{fe}A_{h}$ , the corresponding optimal smoothing factor is

$\mu_{fe,\rm opt}=\min_{\omega}\max_{\boldsymbol{\theta}\in T^{{\rm
high}}}{|1-\omega\widetilde{M}_{fe}\widetilde{A}_{h}|}=\frac{1}{3}\approx
0.333,$ (25)

where the minimum is uniquely achieved at $\omega=\omega_{\rm
opt}=\frac{3}{4}$.

###### Proof.

It can easily be shown that
$\widetilde{M}_{fe}\widetilde{A}_{h}\in\left[\frac{8}{9},\frac{16}{9}\right]$
for $(\theta_{1},\theta_{2})\in T^{\rm high}$. Thus, the optimal $\omega$ is
$\omega=\frac{2}{8/9+16/9}=\frac{3}{4}$. Then, $\mu_{fe,\rm
opt}=1-\frac{2}{16/9+8/9}\frac{8}{9}=\frac{1}{3}$. ∎

From Theorem 3.4, we see that the optimal smoothing factor of 0.333 for the
mass-based relaxation is close to the optimal smoothing factor 0.280 for the
element-wise Vanka patch, and it is better than 0.391 obtained from the
vertex-wise Vanka, see (30), discussed in the next subsection. Thus, mass
matrix could be used as a good approximation to the inverse of the Laplacian
considered here.

####  3.2.2 Vertex-wise Vanka patches in 2D

Now, we analyse the smoothing factor for the vertex-wise patch. Following [6],
the relative Fourier matrix $\Phi$ is

$\Phi=\begin{pmatrix}e^{-\iota\theta_{2}}&0&0&0&0\\\
0&e^{-\iota\theta_{1}}&0&0&0\\\ 0&0&1&0&0\\\ 0&0&0&e^{\iota\theta_{1}}&0\\\
0&0&0&0&e^{\iota\theta_{2}}\end{pmatrix}.$ (26)

The symbol of $M_{v}$ is given by
$\widetilde{M}_{v}=\widetilde{V}^{T}\widetilde{W}\Phi^{T}A_{i}^{-1}\Phi\widetilde{V}$
with

$\widetilde{W}=\frac{1}{5}I,\quad\widetilde{V}=\begin{pmatrix}1\\\ 1\\\ 1\\\
1\\\ 1\end{pmatrix},\quad A_{i}=\frac{1}{h^{2}}\begin{pmatrix}4&0&-1&0&0\\\
0&4&-1&0&0\\\ -1&-1&4&-1&-1\\\ 0&0&-1&4&0\\\ 0&0&-1&0&4\end{pmatrix},$ (27)

where $I$ stands for identity matrix with size $5\times 5$.

It can be shown that

$A_{i}^{-1}=h^{2}\begin{pmatrix}13/48&1/48&1/12&1/48&1/48\\\
1/48&13/48&1/12&1/48&1/48\\\ 1/12&1/12&1/3&1/12&1/12\\\
1/48&1/48&1/12&13/48&1/48\\\ 1/48&1/48&1/12&1/48&13/48\end{pmatrix}.$ (28)

From (26), (27) and (28), we have

$\displaystyle\widetilde{M}_{v}$
$\displaystyle=\widetilde{V}^{T}\widetilde{W}\Phi^{T}A_{i}^{-1}\Phi\widetilde{V},$
$\displaystyle=\frac{h^{2}}{240}\big{(}8(e^{-\iota\theta_{2}}+e^{\iota\theta_{2}}+e^{-\iota\theta_{1}}+e^{\iota\theta_{1}})+68+(e^{-\iota\theta_{2}}+e^{-\iota\theta_{1}})^{2}+(e^{\iota\theta_{2}}+e^{\iota\theta_{1}})^{2}+2(e^{-\iota\theta_{2}}e^{\iota\theta_{1}}+e^{\iota\theta_{2}}e^{-\iota\theta_{1}})\big{)},$
$\displaystyle=\frac{h^{2}}{240}\big{(}16(\cos\theta_{1}+\cos\theta_{2})+68+2\big{(}(\cos\theta_{1}+\cos\theta_{2})^{2}-(\sin\theta_{1}+\sin\theta_{2})^{2}\big{)}+4\cos(\theta_{1}-\theta_{2})\big{)}$
$\displaystyle=\frac{h^{2}}{120}\big{(}(\cos\theta_{1}+\cos\theta_{2}+4)^{2}+(\cos\theta_{1}+\cos\theta_{2})^{2}+16\big{)}.$

Based on the symbol of $M_{v}$, we can write down the corresponding stencil of
$M_{v}$ as follows:

$M_{v}=\frac{h^{2}}{240}\begin{bmatrix}0&0&1&0&0\\\ 0&2&8&2&0\\\ 1&8&68&8&1\\\
0&2&8&2&0\\\ 0&0&1&0&0\end{bmatrix}.$ (29)

Now, we are able to give the optimal smoothing factor for the vertex-wise
Vanka relaxation scheme.

######  Theorem 3.5.

The optimal smoothing factor of $S_{v}$ for the vertex-wise Vanka relaxation
in 2D is

$\mu_{\rm
opt}=\min_{\omega}\max{|1-\omega\widetilde{M}_{v}\widetilde{A}_{h}|}=\frac{9}{23}\approx
0.391,$ (30)

where the minimum is uniquely achieved at $\omega=\omega_{\rm
opt}=\frac{20}{23}\approx 0.8696$.

###### Proof.

We first compute

$\displaystyle\widetilde{M}_{v}\widetilde{A}_{h}$
$\displaystyle=\frac{1}{120}\big{(}(\cos\theta_{1}+\cos\theta_{2}+4)^{2}+(\cos\theta_{1}+\cos\theta_{2})^{2}+16\big{)}\big{(}4-2\cos\theta_{1}-2\cos\theta_{2}\big{)}$
$\displaystyle=\frac{1}{60}\big{(}(4+\eta)^{2}+\eta^{2}+16\big{)}(2-\eta),$

where $\eta=(\cos\theta_{1}+\cos\theta_{2})\in[-2,2]$ with
$(\theta_{1},\theta_{2})\in\left[-\frac{\pi}{2},\frac{3\pi}{2}\right]^{2}$.

Let $g(\eta)=\frac{1}{60}\Big{(}(4+\eta)^{2}+\eta^{2}+16\Big{)}(2-\eta)$. We
find that

$g^{\prime}(\eta)=\frac{1}{30}(-3\eta^{2}-4\eta-8)<0,\quad\forall\eta\in\mathbb{R}.$

This means that $g(\eta)$ is a decreasing function. Thus, for $\eta\in[-2,2]$,
we have

$0=g(\eta=2)\leq g(\eta)\leq g(\eta=-2)=\frac{8}{5}.$

This means that $g(\eta)\in[0,\frac{8}{5}]$.

If we restrict $(\theta_{1},\theta_{2})\in T^{\rm high}$, then
$\eta\in[-2,1]$. In this situation, we have

$\frac{7}{10}=g(\eta=1)\leq g(\eta)\leq g(\eta=-2)=\frac{8}{5}.$

Since $\widetilde{M}_{v}\widetilde{A}_{h}\in[\frac{7}{10},\frac{8}{5}]$ for
$(\theta_{1},\theta_{2})\in T^{\rm high}$, we have $\omega_{\rm
opt}=\frac{2}{7/10+8/5}=\frac{20}{23}$ and $\mu_{\rm
opt}=1-\frac{20}{23}\frac{7}{10}=\frac{9}{23}$. ∎

######  Remark 3.1.

Note that in 2D the element-wise Vanka, see (21), uses fewer points than that
of vertex-wise Vanka, see (29). However, the corresponding optimal smoothing
factor of element-wise Vanka, see (19), is smaller than that of vertex-wise
Vanka, see (30), which is different than the case in 1D.

###  3.3 Extension to the 3D case

While we do not include a smoothing analysis of Vanka-type solvers for the 3D
case, we can still make a few interesting observations. In particular,
motivated by our findings on the potential role of the mass matrix for
relaxation, we further explore the scaled mass matrix in 3D as an
approximation to the inverse of the Laplacian. Let

$M_{3}=h^{-4}M_{e}\otimes M_{e}\otimes M_{e},$

where $M_{e}$ is defined in (11). The symbol of $M_{3}$ can be obtained by
tensor product given by

$\widetilde{M}_{3}=\frac{h^{2}}{27}(2+\cos\theta_{1})(2+\cos\theta_{2})(2+\cos\theta_{3}).$

The symbol of $A_{h}$ in 3D is

$\widetilde{A}_{h}=\frac{1}{h^{2}}2(3-\cos\theta_{1}-\cos\theta_{2}-\cos\theta_{3}).$

Let $M_{J}=\frac{6}{h^{2}}I$, which is the Jacobi matrix.

######  Theorem 3.6.

If we consider the point-wise damped Jacobi as a preconditioner for the
Laplacian, then the corresponding optimal smoothing factor of $S=I-\omega
M_{J}A_{h}$ is

$\mu_{J,\rm opt}=\min_{\omega}\max_{\boldsymbol{\theta}\in T^{{\rm
high}}}{|1-\omega\widetilde{M}_{J}\widetilde{A}_{h}|}=\frac{5}{7}\approx
0.714,$ (31)

where the minimum is uniquely achieved at $\omega_{\rm opt}=\frac{6}{7}\approx
0.857$.

###### Proof.

Note that
$\widetilde{M}_{J}\widetilde{A}_{h}=\frac{1}{3}(3-\cos\theta_{1}-\cos\theta_{2}-\cos\theta_{3})$.
For $\boldsymbol{\theta}\in T^{\rm high}$,
$\widetilde{M}_{J}\widetilde{A}_{h}\in[1/3,2]$. Thus, to minimize
$\mu(\omega)=|1-\omega\widetilde{M}_{J}\widetilde{A}_{h}|$, we require
$\omega=\frac{2}{1/3+2}=\frac{6}{7}$. It follows $\mu_{J,\rm
opt}=1-\frac{6}{7}\frac{1}{3}=\frac{5}{7}$. ∎

Next, we consider mass-based relaxation scheme for the Laplacian in 3D.

######  Theorem 3.7.

Given the scaled mass stencil $M_{3}$ in (3.3) and mass-based relaxation
scheme $S_{m}=I-\omega M_{3}A_{h}$ in 3D, the corresponding optimal smoothing
factor is

$\mu_{m,\rm opt}=\min_{\omega}\max_{\boldsymbol{\theta}\in T^{{\rm
high}}}{|1-\omega\widetilde{M}_{3}\widetilde{A}_{h}|}=\frac{131}{212}\approx
0.618,$ (32)

where the minimum is uniquely achieved at $\omega=\omega_{\rm
opt}=\frac{729}{848}\approx 0.860$.

###### Proof.

Let

$T=\widetilde{M}_{3}\widetilde{A}_{h}=\frac{2}{27}(2+\cos\theta_{1})(2+\cos\theta_{2})(2+\cos\theta_{2})(3-\cos\theta_{1}-\cos\theta_{2}-\cos\theta_{3}).$

Define $g(x,y,z)=\frac{2}{27}(2+x)(2+y)(2+z)(3-x-y-z)$ with $x,y,z\in[-1,1]$.
To find the extremal values of $g$, we will consider its derivatives and the
function values at the boundary of underlying domain. We compute the
derivatives of $g$ with respect to $x,y$ and $z$, given by

$\displaystyle g_{x}$ $\displaystyle=\frac{2}{27}(2+y)(2+z)(1-2x-y-z),$
$\displaystyle g_{y}$ $\displaystyle=\frac{2}{27}(2+x)(2+z)(1-2y-x-z),$
$\displaystyle g_{z}$ $\displaystyle=\frac{2}{27}(2+x)(2+y)(1-2z-x-y).$

Solving $g_{x}=g_{y}=g_{z}=0$ with $x,y,z\in[-1,1]$ gives $x=y=z=1/4$.
However, $\boldsymbol{\theta^{*}}=(\theta_{1},\theta_{2},\theta_{3})$ such
that $(\cos\theta_{1},\cos\theta_{2},\cos\theta_{3})=(1/4,1/4,1/4)$ does not
belong to $T^{\rm high}$.

Let us define $\Omega_{1}=[-1,1]^{3}$, $\Omega_{2}=[0,1]^{3}$, and
$\Omega=\Omega_{1}\big{\backslash}\Omega_{2}$. Note that $\Omega$ corresponds
to $\boldsymbol{\theta}\in T^{\rm high}$. To find the extremal values of $T$
for $\boldsymbol{\theta}\in T^{\rm high}$, we only need to find the extremal
values of $T$ at the boundary of $\Omega$, denoted as $\partial\Omega$. Note
that $\partial\Omega$ contains the following four cases.
Case 1:

$\displaystyle x=-1,(y,z)\in[-1,1]^{2},$ $\displaystyle
y=-1,(x,z)\in[-1,1]^{2},$ $\displaystyle z=-1,(x,y)\in[-1,1]^{2}.$

Case 2:

$\displaystyle x=1,(y,z)\in[-1,0]\times[-1,1],$ $\displaystyle
y=1,(x,z)\in[-1,0]\times[-1,1],$ $\displaystyle
z=1,(x,y)\in[-1,0]\times[-1,1].$

Case 3:

$\displaystyle x=0,(y,z)\in[0,1]^{2},$ $\displaystyle y=0,(x,z)\in[0,1]^{2}$
$\displaystyle z=0,(x,y)\in[0,1]^{2}.$

Case 4:

$\displaystyle x=1,(y,z)\in[0,1]\times[-1,0],$ $\displaystyle
y=1,(x,z)\in[0,1]\times[-1,0],$ $\displaystyle z=1,(x,y)\in[0,1]\times[-1,0].$

Due to the symmetry of $g(x,y,z)$ and our interest of maximum and minimum of
$g(x,y,z)$, we only need to consider the following sets

$\displaystyle\mathcal{D}_{1}$
$\displaystyle=\left\\{(x,y,z)|\,x=-1,(y,z)\in[-1,1]^{2}\right\\},$
$\displaystyle\mathcal{D}_{2}$
$\displaystyle=\left\\{(x,y,z)|\,x=1,(y,z)\in[-1,0]\times[-1,1]\right\\},$
$\displaystyle\mathcal{D}_{3}$
$\displaystyle=\left\\{(x,y,z)|\,x=0,(y,z)\in[0,1]^{2}\right\\},$
$\displaystyle\mathcal{D}_{4}$
$\displaystyle=\left\\{(x,y,z)|\,x=1,(y,z)\in[0,1]\times[-1,0]\right\\}.$

For $\boldsymbol{\theta}\in T^{\rm high}$, we check the extremal values of
$g(x,y,z)$ on the sets
$\mathcal{D}_{1}\bigcup\mathcal{D}_{2}\bigcup\mathcal{D}_{3}\bigcup\mathcal{D}_{4}$,
and find that the maximum of $T$ is $\frac{4\cdot 7^{3}}{3^{6}}$ achieved at
$(\cos\theta_{1},\cos\theta_{2},\cos\theta_{3})=(0,1/3,1/3)$ and the smallest
value of $T$ is $\frac{4}{9}$ with
$(\cos\theta_{1},\cos\theta_{2},\cos\theta_{3})=(-1,-1,-1)$. Thus, the optimal
parameter is $\omega=\frac{2}{4/9+4\cdot 7^{3}/3^{6}}=\frac{729}{848}\approx
0.860$ and the corresponding smoothing factor is $\mu_{m,\rm
opt}=1-\frac{729}{848}\frac{4}{9}=\frac{131}{212}\approx 0.618$. ∎

######  Remark 3.2.

One might consider a tensor-product generalization, i.e.,
$M_{v,3}=h^{-4}M_{v}\otimes M_{v}\otimes M_{v}$, where $M_{v}$ is defined in
(15), as an approximation to the inverse of the Laplacian in 3D. However, with
this choice, we find that the numerically optimal smoothing factor is 0.800,
which is larger than 0.618, see (32), obtained from the mass-based relaxation.
Thus, we do not further explore $M_{v,3}$.

##  4 Numerical experiments

In this section we compute the smoothing factor to validate our theoretical
results. We also report the LFA two-grid convergence factor and compare it to
the corresponding smoothing factor.

In general, the two-grid error-propagation operator can be expressed as [14,
18]

$E=S^{\nu_{2}}(I-PA_{H}^{-1}RA_{h})S^{\nu_{1}},$ (33)

where $R$ is the restriction operator from grid $h$ to grid $H$, $P$ is the
interpolation operator from grid $H$ to grid $h$, $A_{H}$ represents the
coarse-grid operator, and the integers $\nu_{1}$ and $\nu_{2}$ are the numbers
of pre- and post-relaxation sweeps, respectively. Here, $S$ is the relaxation
error operator defined in (6).

######  Definition 4.1.

The LFA two-grid convergence factor [15, 19] for $E$ is defined as

$\rho=\max_{\boldsymbol{\theta}\in T^{\rm
low}}\left\\{\rho(\widetilde{E}(\boldsymbol{\theta},\omega))\,\right\\},$ (34)

where $\omega$ is a parameter, $\widetilde{E}$ denotes the symbol of the two-
grid error operator $E$, and $\rho(\widetilde{E}(\boldsymbol{\theta},\omega))$
denotes the spectral radius of the matrix $\widetilde{E}$.

For more details on how to obtain the symbol of the two-grid error-propagation
operator $E$, see [15, 19]. In our tests, we consider $P$ to be the standard
(bi)linear interpolation and take $R=P^{T}$ and $A_{H}=RA_{h}P$. We take
$\rho$ to be the LFA prediction sampled at 64 equispaced points in each
dimension of the Fourier domain. For simplicity, we denote by $T(\nu)$ the
two-grid method with $\nu=\nu_{1}+\nu_{2}$.

In Figure 2 we present the eigenvalue distribution of $\widetilde{E}$ for the
Vanka-type method in 2D with the optimal value of $\omega$. We see that all
eigenvalues of $\widetilde{E}$ are real and their largest magnitude matches
the optimal smoothing factor.

Figure 3 shows the magnitude of eigenvalues of $\widetilde{E}_{e}$ and
$\widetilde{S}_{v}$ as a function of the Fourier modes,
$(\theta_{1},\theta_{2})$. We see that the eigenvalues of the two-grid error
operator are distributed almost evenly in the entire Fourier domain. On the
other hand, for the vertex-wise Vanka, the largest magnitude occurs at the
$(0,\pm\pi)$ or $(\pm\pi,0)$, see Figure 4. From the right of Figures 3 and 4,
we see that the relaxation scheme reduces the high frequency errors.

Figure 2: Eigenvalue distribution of $\widetilde{E}$ for Vanka-type method
with $\nu=1$. Left: element-wise Vanka. Right: vertex-wise Vanka.

Figure 3:  Element-wise Vanka with $\nu=1$. Left: the magnitude of eigenvalues
of $\widetilde{E}_{e}$ as a function of the Fourier modes,
$(\theta_{1},\theta_{2})$. Right: the magnitude of eigenvalues of
$\widetilde{S}_{e}$ as a function of the Fourier modes,
$(\theta_{1},\theta_{2})$.

Figure 4:  Vertex-wise Vanka with $\nu=1$. Left: the magnitude of eigenvalues of $\widetilde{E}_{v}$ as a function of the Fourier modes, $(\theta_{1},\theta_{2})$. Right: the magnitude of eigenvalues of $\widetilde{S}_{v}$ as a function of the Fourier modes, $(\theta_{1},\theta_{2})$. Table 2: LFA predicted two-grid convergence factor, ${\rho}$ with $\omega_{\rm opt}$ obtained from our analytical result and $h=\frac{1}{64}$. Cycle | $\omega_{\rm opt}$ | $\mu_{\rm opt}$ | $TG(1)$ | $TG(2)$ | $TG(3)$ | $TG(4)$
---|---|---|---|---|---|---
1D
$\rho_{e}$ | 12/17 | 0.059 | 0.059 | 0.059 | 0.040 | 0.031
$\rho_{v}$ | 81/104 | 0.038 | 0.091 | 0.033 | 0.022 | 0.017
2D
$\rho_{e}$ | 24/25 | 0.280 | 0.280 | 0.092 | 0.059 | 0.045
$\rho_{v}$ | 20/23 | 0.391 | 0.391 | 0.153 | 0.076 | 0.055
$\rho_{fe}$ | 3/4 | 0.333 | 0.333 | 0.111 | 0.037 | 0.029
3D
$\rho_{J}$ | 6/7 | 0.714 | 0.714 | 0.510 | 0.364 | 0.260
$\rho_{m}$ | 729/848 | 0.618 | 0.618 | 0.382 | 0.236 | 0.146

From Table 2, we see that the smoothing factors match the two-grid convergence
factor with $\nu=1$ (except for the vertex-wise patch for 1D), which is as
expected since the smoothing factor typically offers a sharp prediction of the
two-grid convergence factor. As $\nu$ increases, we see a little degradation
of $\rho(\nu)$, that is $\rho(\nu)>\mu^{\nu}$ for Vanka-type relaxation, and
so does the mass-based relaxation scheme with $\nu=4$.

To explore the gap between the smoothing factor and two-grid convergence
factor for the vertex-wise Vanka in 1D, we run our LFA code to minimize the
LFA two-grid convergence factor with respect to parameter $\omega$ by an
exhaustive search with stepsize 0.02, and find that the optimal LFA two-grid
convergence factor is 0.067 with $\omega=0.80$ and the corresponding smoothing
factor is the same as the two-grid LFA convergence factor. For the choice of
$\omega=\frac{81}{104}$ that minimizes the smoothing factor, the gap between
the LFA two-grid convergence factor and the smoothing factor is reasonable, as
is noted in the literature; see, for example [8]. The smoothing factor is a
sharp prediction of the two-grid convergence factor.

######  Remark 4.1.

The LFA two-grid convergence and smoothing factors seem independent of the
meshsize $h$. We have tested different values of $h$ for Table 2 to confirm
this. Further details are omitted.

##  5 Conclusions

We have presented a theoretical analysis of the optimal multigrid smoothing
factor for two types of additive Vanka smoothers, applied to the Poisson
equation discretized by the standard centered finite difference scheme. The
smoothers are shown to be highly efficient. We have found that the element-
wise Vanka is closely related to the mass matrix obtained from the (bi)linear
finite element method. This observation propels us to use the mass matrix as
an approximation to the inverse of the Laplacian, and this yields rapid
convergence. It is an interesting connection between the finite element method
and finite difference methods.

Solvers for equations related to the Laplacian are a natural first step in
developing new algorithms for more complex problems, such as the Stokes
equations, Navier-Stokes equations and other saddle-point problems. Using a
similar approach as part of the development of fast solvers for those problems
may prove computationally beneficial.

## References

  * [1] J. H. Adler, T. R. Benson, and S. P. MacLachlan.  Preconditioning a mass-conserving discontinuous Galerkin discretization of the Stokes equations.  Numerical Linear Algebra with Applications, 24(3):e2047, 2017.
  * [2] J. H. Adler, Y. He, X. Hu, S. MacLachlan, and P. Ohm.  Monolithic multigrid for a reduced-quadrature discretization of poroelasticity.  arXiv preprint arXiv:2107.04060, 2021.
  * [3] L. Claus and M. Bolten.  Nonoverlapping block smoothers for the Stokes equations.  Numerical Linear Algebra with Applications, page e2389, 2021.
  * [4] A. P. de la Riva, C. Rodrigo, and F. J. Gaspar.  A Robust Multigrid Solver for Isogeometric Analysis Based on Multiplicative Schwarz Smoothers.  SIAM J. Sci. Comput., 41(5):S321–S345, 2019.
  * [5] Á. P. de la Riva, C. Rodrigo, and F. J. Gaspar.  A two-level method for isogeometric discretizations based on multiplicative Schwarz iterations.  Computers & Mathematics with Applications, 100:41–50, 2021.
  * [6] P. E. Farrell, Y. He, and S. P. MacLachlan.  A local fourier analysis of additive Vanka relaxation for the Stokes equations.  Numerical Linear Algebra with Applications, 28(3):e2306, 2021.
  * [7] S. Franco, C. Rodrigo, F. Gaspar, and M. Pinto.  A multigrid waveform relaxation method for solving the poroelasticity equations.  Computational and Applied Mathematics, pages 1–16, 2018.
  * [8] Y. He and S. MacLachlan.  Two-level Fourier analysis of multigrid for higher-order finite-element discretizations of the Laplacian.  Numerical Linear Algebra with Applications, 27(3):e2285, 2020.
  * [9] Y. He, S. Rhebergen, and H. De Sterck.  Local Fourier analysis of multigrid for hybridized and embedded discontinuous Galerkin methods.  SIAM Journal on Scientific Computing, (0):S612–S636, 2021.
  * [10] V. John.  Higher order finite element methods and multigrid solvers in a benchmark problem for the 3D Navier–Stokes equations.  International Journal for Numerical Methods in Fluids, 40(6):775–798, 2002.
  * [11] V. John and L. Tobiska.  Numerical performance of smoothers in coupled multigrid methods for the parallel solution of the incompressible Navier-Stokes equations.  International Journal For Numerical Methods In Fluids, 33(4):453–473, Jan 2000.
  * [12] S. Manservisi.  Numerical analysis of Vanka-type solvers for steady Stokes and Navier-Stokes flows.  SIAM Journal on Numerical Analysis, 44(5):2025–2056, 2006.
  * [13] S. Saberi, G. Meschke, and A. Vogel.  A restricted additive Vanka smoother for geometric multigrid.  arXiv preprint arXiv:2103.10127, 2021.
  * [14] K. Stüben and U. Trottenberg.  Multigrid methods: Fundamental algorithms, model problem analysis and applications.  In Multigrid methods, pages 1–176. Springer, 1982.
  * [15] U. Trottenberg, C. W. Oosterlee, and A. Schüller.  Multigrid.  Academic Press, Inc., San Diego, CA, 2001.
  * [16] S. P. Vanka.  Block-implicit multigrid solution of Navier-Stokes equations in primitive variables.  Journal of Computational Physics, 65:138–158, 1986.
  * [17] A. Voronin, Y. He, S. MacLachlan, L. N. Olson, and R. Tuminaro.  Low-order preconditioning of the Stokes equations.  Numerical Linear Algebra with Applications, 2021.  To appear.
  * [18] P. Wesseling.  An introduction to multigrid methods.  Pure and Applied Mathematics (New York). John Wiley & Sons, Ltd., Chichester, 1992.
  * [19] R. Wienands and W. Joppich.  Practical Fourier analysis for multigrid methods.  CRC press, 2004.