#  _What Are You Token About?_
Dense Retrieval as Distributions Over the Vocabulary

Ori Ram1 Liat Bezalel1 Adi Zicher1
Yonatan Belinkov2 Jonathan Berant1 Amir Globerson1
1Blavatnik School of Computer Science, Tel Aviv University 2Technion – IIT,
Israel
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Supported by the Viterbi Fellowship in the Center for Computer Engineering at
the Technion.

###### Abstract

Dual encoders are now the dominant architecture for dense retrieval. Yet, we
have little understanding of how they represent text, and why this leads to
good performance. In this work, we shed light on this question via
distributions over the vocabulary. We propose to interpret the vector
representations produced by dual encoders by projecting them into the model’s
vocabulary space. We show that the resulting projections contain rich semantic
information, and draw connection between them and sparse retrieval. We find
that this view can offer an explanation for some of the failure cases of dense
retrievers. For example, we observe that the inability of models to handle
tail entities is correlated with a tendency of the token distributions to
forget some of the tokens of those entities. We leverage this insight and
propose a simple way to enrich query and passage representations with lexical
information at _inference_ time, and show that this significantly improves
performance compared to the original model in zero-shot settings, and
specifically on the BEIR benchmark.111Our code is publicly available at
https://github.com/oriram/dense-retrieval-projections.

Figure 1:  An example of our framework. We run the question “Where was Michael
Jack born?” through the question encoder of DPR Karpukhin et al. (2020), and
project the question representation $\bm{e}_{q}$ to the vocabulary space using
BERT’s masked language modeling head Devlin et al. (2019). The result is a
distribution over the vocabulary, $Q$. We apply the same procedure for
passages as well. These projections enable reasoning about and improving
retrieval representations.  Figure 2: A success case from Natural Questions
(top) and a failure case from EntityQuestions (bottom) of DPR Karpukhin et al.
(2020), explained via projecting question and (its relevant) passage
representations to the vocabulary space. Tokens in the top-20 of both question
and passage vocabulary projections are marked in bold.

##  1 Introduction

Dense retrieval models based on neural text representations have proven very
effective Karpukhin et al. (2020); Qu et al. (2021); Ram et al. (2022);
Izacard et al. (2022a, b), improving upon strong traditional sparse models
like BM25 Robertson and Zaragoza (2009). However, when applied off-the-shelf
(i.e., in out-of-domain settings) they often experience a severe drop in
performance Thakur et al. (2021); Sciavolino et al. (2021); Reddy et al.
(2021). Moreover, the reasons for such failures are poorly understood, as the
information captured in their representations remains under-investigated.

In this work, we present a new approach for interpreting and reasoning about
dense retrievers, through distributions induced by their query222Throughout
the paper, we use query and question interchangeably. and passage
representations when projected to the vocabulary space, namely distributions
over their vocabulary space (Figure 1). Such distributions enable a better
understanding of the representational nature of dense models and their
failures, which paves the way to simple solutions that improve their
performance.

We begin by showing that dense retrieval representations can be projected to
the vocabulary space, by feeding them through the masked language modeling
(MLM) head of the pretrained model they were initialized from without any
further training. This operation results in distributions over the vocabulary,
which we refer to as query vocabulary projections and passage vocabulary
projections.

Surprisingly, we find these projections to be highly interpretable to humans
(Figure 2; Table 1). We analyze these projections and draw interesting
connections between them and well-known concepts from sparse retrieval (§5).
First, we highlight the high coverage of tokens shared by the query and the
passage in the top-$k$ of their projections. This obersvation suggests that
the lexical overlap between query and passages plays an important role in the
retrieval mechanism. Second, we show that vocabulary projections of passages
they are likely to contain words that appear in queries about the given
passage. Thus, they can be viewed as predicting the questions one would ask
about the passage. Last, we show that the model implicitly implements query
expansion Rocchio (1971). For example, in Figure 2 the query is “How many
judges currently serve on the Supreme court?”, and the words in the query
projection $Q$ include “justices” (the common way to refer to them) and “nine”
(the correct answer).

The above findings are especially surprising due to the fact that these
retrieval models are fine-tuned in a contrastive fashion, and thus do not
perform any prediction over the vocabulary or make any use of their language
modeling head during fine-tuning. In addition, these representations are the
result of running a deep transformer network that can implement highly complex
functions. Nonetheless, model outputs remain “faithful” to the original
lexical space learned during pretraining.

We further show that our approach is able to shed light on the reasons for
which dense retrievers struggle with simple entity-centric questions
Sciavolino et al. (2021). Through the lens of vocabulary projections, we
identify an interesting phenomenon: dense retrievers tend to “ignore” some of
the tokens appearing in a given passage. This is reflected in the ranking
assigned to such tokens in the passage projection. For example, the word
“michael” in the bottom example of Figure 2 is ranked relatively low (even
though it appears in the passage title), thereby hindering the model from
retrieving this passage. We refer to this syndrome as token amnesia (§6).

We leverage this insight and suggest a simple inference-time fix that enriches
dense representations with lexical information, addressing token amnesia. We
show that lexical enrichment significantly improves performance compared to
vanilla models on the challenging BEIR benchmark Thakur et al. (2021) and
additional datasets. For example, we boost the performance of the strong MPNet
model on BEIR from 43.1% to 44.1%.

Taken together, our analyses and results demonstrate the great potential of
vocabulary projections as a framework for more principled research and
development of dense retrieval models.

##  2 Background

In this work, we suggest a simple framework for interpreting dense retrieves,
via projecting their representations to the vocabulary space. This is done
using the (masked) language modeling head of their corresponding pretrained
model. We begin by providing the relevant background.

###  2.1 Masked Language Modeling

Most language models based on encoder-only transformers Vaswani et al. (2017)
are pretrained using some variant of the masked language modeling (MLM) task
Devlin et al. (2019); Liu et al. (2019); Song et al. (2020), which involves
masking some input tokens, and letting the model reconstruct them.

Specifically, for an input sequence $x_{1},...,x_{n}$, the transformer encoder
is applied to output contextualized token representations
$\bm{h}_{1},...,\bm{h}_{n}\in\mathbb{R}^{d}$. Then, to predict the missing
tokens, an MLM head is applied to their contextualized representations. The
MLM head is a function that takes a vector $\bm{h}\in\mathbb{R}^{d}$ as input
and returns a distribution $P$ over the model’s vocabulary $\mathcal{V}$,
defined as follows:

$\begin{split}\text{MLM-
Head}(\bm{h})[i]=\frac{\exp(\bm{v}_{i}^{\top}g(\bm{h}))}{\sum_{j\in\mathcal{V}}\exp(\bm{v}_{j}^{\top}g(\bm{h}))}\end{split}$
(1)

$g:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ is a potentially non-linear
function (e.g., a fully connected layer followed by a LayerNorm for BERT;
Devlin et al. 2019), and $\bm{v}_{i}\in\mathbb{R}^{d}$ corresponds to the
static embedding of the $i$-th item in the vocabulary.

###  2.2 Dense Retrieval

In dense retrieval, we are given a corpus of passages
$\mathcal{C}=\\{p_{1},...,p_{m}\\}$ and a query $q$ (e.g., a question or a
fact to check), and we wish to compute query and passage representations
($\bm{e}_{q}$ and $\bm{e}_{p}$, respectively) such that similarity in this
space implies high relevance of a passage to the query. Formally, let
$\text{Enc}_{Q}$ be a query encoder and $\text{Enc}_{P}$ a passage encoder.
These encoders are mappings from the input text to a vector in
$\mathbb{R}^{d}$, and are obtained by fine-tuning a given LLM. Specifically,
they return a pooled version of the LLM contextualized embeddings (e.g., the
[CLS] embedding or mean pooling). We denote the embedding of the query and
passage vectors as follows:

$\begin{split}\bm{e}_{q}&=\text{Enc}_{Q}(q)\\\
\bm{e}_{p}&=\text{Enc}_{P}(p)\\\ \end{split}$ (2)

To fine-tune retrievers, a similarity measure $s(q,p)$ is defined (e.g., the
dot-product between $\bm{e}_{q}$ and $\bm{e}_{q}$ or their cosine similarity)
and the model is trained in a contrastive manner to maximize retriever
accuracy Lee et al. (2019); Karpukhin et al. (2020). Importantly, in this
process, the MLM head function does not change at all.

##  3 Vocabulary Projections

We now describe our framework for projecting query and passage representations
of dense retrievers to the vocabulary space. Given a dense retrieval model, we
utilize the MLM head of the model it was initialized from to map from encoder
output representations to distributions over the vocabulary (Eq. 1). For
example, for DPR Karpukhin et al. (2020) we take BERT’s MLM head, as DPR was
initialized from BERT. Given a query $q$, we use the query encoder
$\text{Enc}_{Q}$ to obtain its representation $\bm{e}_{q}$ as in Eq. 2.
Similarly, for a passage $p$ we apply the passage encoder $\text{Enc}_{P}$ to
get $\bm{e}_{p}$. We then apply the MLM head as in Eq. (1) to obtain the
vocabulary projection:

$\begin{split}Q&=\text{MLM-Head}(\bm{e}_{q})\\\ P&=\text{MLM-
Head}(\bm{e}_{p})\end{split}$ (3)

Note that it is not clear a-priori that $Q$ and $P$ will be meaningful in any
way, as the encoder model has been changed since pretraining, while the MLM-
head function remains fixed. Moreover, the MLM function has not been trained
to decode “pooled” sequence-level representations (i.e., the results of CLS or
mean pooling) during pretraining. Despite this intuition, in this work we
argue that $P$ and $Q$ are actually highly intuitive and can facilitate a
better understanding of dense retrievers.

Question |  top-20 in $Q$ |  Passage |  top-20 in $P$
---|---|---|---
where do the great lakes meet the ocean (A: the saint lawrence river) |  lakes lake shore ocean confluence river water north canada meet east land rivers canoe sea border michigan connecting both shores |  the great lakes , also called the laurent ##ian great lakes and the great lakes of north america , are a series of inter ##connected freshwater lakes located primarily in the upper mid - east region of north america , on the canada – united states border , which connect to the atlantic ocean through the saint lawrence river . they consist of lakes superior , michigan , huron … |  lakes lake the canada great freshwater water region ontario these central river rivers large basin core area erie all four
southern soul was considered the sound of what independent record label (A: motown) |  southern music label soul motown blues nashville vinyl sound independent labels country records genre dixie record released gospel jazz south |  soul music . the key sub ##gen ##res of soul include the detroit ( motown ) style , a rhythmic music influenced by gospel ; " deep soul " and " southern soul " , driving , energetic soul styles combining r & b with southern gospel music sound ; … which came out of the rhythm and blues style … |  soul music jazz funk blues rock musical fusion genre black pure classical genres pop southern melody art like rich urban
who sings does he love me with re ##ba (A: linda davis) |  duet song love music solo re he motown me his " pa album songs honey reprise bobby i peggy blues |  " does he love you " is a song written by sandy knox and billy st ##rit ##ch , and recorded as a duet by american country music artists re ##ba mc ##ent ##ire and linda davis … |  he you him i it she his john we love paul who me does did yes why they how this
Table 1: Examples of questions and gold passages from the development set of
Natural Questions, along with their 20 top-scored tokens in projections of DPR
representations. Green tokens represent the lexical overlap signal (i.e.,
tokens that appear in both the question and the passage). Blue tokens
represent query expansion (i.e., tokens that do not appear in the question but
do appear in the passage).

##  4 Experiment Setup

To evaluate our framework and method quantitatively, we consider several dense
retrieval models and datasets.

###  4.1 Models

We now list the retrievers used to demonstrate our framework and method. All
dense models share the same architecture and size (i.e., that of BERT-base;
110M parameters), and all were trained in a contrastive fashion with in-batch
negatives—the prominent paradigm for training dense models Lee et al. (2019);
Karpukhin et al. (2020); Chang et al. (2020); Qu et al. (2021); Ram et al.
(2022); Izacard et al. (2022a); Ni et al. (2022); Chen et al. (2022). For the
analysis, we use DPR Karpukhin et al. (2020) and BERT Devlin et al. (2019) as
its pretrained baseline. For the results of our method, we also use S-MPNet
Reimers and Gurevych (2019) and Spider Ram et al. (2022). Our sparse retrieval
model is BM25 Robertson and Zaragoza (2009). We refer the reader to App. A for
more details.

###  4.2 Datasets

We follow prior work Karpukhin et al. (2020); Ram et al. (2022) and consider
six common open-domain question answering (QA) datasets for the evaluation of
our framework: Natural Questions (NQ; Kwiatkowski et al. 2019), TriviaQA Joshi
et al. (2017), WebQuestions (WQ; Berant et al. 2013), CuratedTREC (TREC;
Baudiš and Šedivý 2015), SQuAD Rajpurkar et al. (2016) and EntityQuestions
(EntityQs; Sciavolino et al. 2021). We also consider the BEIR Thakur et al.
(2021) and the MTEB Muennighoff et al. (2022) benchmarks.

###  4.3 Implementation Details

Our code is based on the official repository of DPR Karpukhin et al. (2020),
built on Hugging Face Transformers Wolf et al. (2020).

For the six QA datasets, we use the Wikipedia corpus standardized by Karpukhin
et al. (2020), which contains roughly 21 million passages of a hundred words
each. For dense retrieval over this corpus, we apply exact search using FAISS
Johnson et al. (2021). For sparse retrieval we use Pyserini Lin et al. (2021).

##  5 Analyzing Dense Retrievers via Vocabulary Projections

In Section 3, we introduce a new framework for interpreting representations
produced by dense retrievers. Next, we describe empirical findings that shed
new light on what is encoded in these representations. Via vocabulary
projections, we draw connections between dense retrieval and well-known
concepts from sparse retrieval like lexical overlap (§5.1), query prediction
(§5.2) and query expansion (§5.3).

###  5.1 The Dominance of Lexical Overlap

Figure 3:  The percentage of tokens shared by questions (from the development
set of NQ) and their gold passages (i.e., the lexical overlap signal) that are
covered by the top-$k$ tokens of question vocabulary projection $Q$ and
passage vocabulary projection $P$ as a function of $k$. Stop words and
punctuation marks are excluded from this analysis.

Tokens shared by questions and their corresponding gold passages constitute
the lexical overlap signal in retrieval, used by sparse models like BM25. We
start by asking: how prominent are they in vocabulary projections? Figure 3
illustrates the coverage of these tokens in $Q$ and $P$ for DPR after
training, compared to its initialization before training (i.e., BERT with mean
or CLS pooling). In other words, for each $k$ we check what is the percentage
of shared tokens ranked in the top-$k$ of $Q$ and $P$. Results suggest that
after training, the model learns to rank shared tokens much higher than
before. Concretely, 63% and 53% of the shared tokens appear in the top-20
tokens of $Q$ and $P$ respectively, compared to only 16% and 8% in BERT (i.e.,
before training). These numbers increase to 78% and 69% of the shared tokens
that appear in the top-100 tokens of $Q$ and $P$. In addition, we observed
that for 71% of the questions, the top-scored token in $Q$ appears in both the
question and the passage (App. B). These findings suggest that even for dense
retrievers—which do not operate at the lexical level—lexical overlap remains a
highly dominant signal.

###  5.2 Passage Encoders as Query Prediction

Our next analysis concerns the role of passage encoders. In §5.1, we show that
tokens shared by the question and its gold passage are ranked high in both $Q$
and $P$. However, passages contain many tokens, and the shared tokens
constitute only a small fraction of all tokens. We hypothesize that out of
passage tokens, those that are likely to appear in relevant questions receive
higher scores in $P$ than others. If this indeed the case, it implies that
passage encoders implicitly learn to predict which of the passage tokens will
appear in relevant questions. To test our hypothesis, we analyze the ranks of
question and passage tokens in passage vocabulary projections, $P$. Formally,
let $\mathcal{T}_{q}$ and $\mathcal{T}_{p}$ be the sets of tokens in a
question $q$ and its gold passage $p$, respectively. Table 2 shows the token-
level mean reciprocal rank (MRR) of these sets in $P$. We observe that tokens
shared by $q$ and $p$ (i.e., $\mathcal{T}_{q}\cap\mathcal{T}_{p}$) are ranked
significantly higher than other passage tokens (i.e., $\mathcal{T}_{p}$). For
example, in DPR the MRR of shared tokens is 26.1, while that of other passage
tokens is only 3.0. In addition, the MRR of shared tokens in BERT is only 1.4.
These findings support our claim that tokens that appear in relevant questions
are ranked higher than others, and that this behavior is acquired during fine-
tuning.

|  | Token-Level MRR in $P$
---|---|---
|  |  DPR | BERT (mean)
Passage tokens | $\mathcal{T}_{p}$ |  3.0 |  0.5
Question tokens | $\mathcal{T}_{q}$ | 17.3 |  1.0
Shared tokens | $\mathcal{T}_{q}\cap\mathcal{T}_{p}$ | 26.1 |  1.4
Table 2: An analysis of token-level MRR (in %) in passage vocabulary
projections $P$ on the development set of NQ. For a question $q$ and its gold
positive passage $p$, $\mathcal{T}_{q}$ and $\mathcal{T}_{p}$ are the
corresponding sets of tokens, excluding stop words and punctuations. For a set
$\mathcal{T}$, we report
$\frac{1}{|\mathcal{T}|}\sum_{t\in\mathcal{T}}\frac{1}{\text{rank}_{P}(t)}$.

###  5.3 Query Encoders Implement Query Expansion

To overcome the “vocabulary mismatch” problem (i.e., when question-document
pairs are semantically relevant, but lack significant lexical overlap), query
expansion methods have been studied extensively Rocchio (1971); Voorhees
(1994); Zhao and Callan (2012); Mao et al. (2021). The main idea is to expand
the query with additional terms that will better guide the retrieval process.
We define a token as a query expansion if it does not appear in the query
itself but does appear in the query projection $Q$, and also in the gold
passage of that query $p$ (excluding stop words and punctuation marks). Figure
4 shows the percentage of queries with at least one query expansion token in
the top-$k$ as a function of $k$ for DPR and the BERT baseline (i.e., before
DPR training). We observe that after training, the model promotes query
expansion tokens to higher ranks than before. In addition, we found that
almost 14% of the tokens in the top-5 of $Q$ are query expansion tokens (cf.
App B).

We note that there are two interesting classes of query expansion tokens: (1)
synonyms of question tokens, as well as tokens that share similar semantics
with tokens in $q$ (e.g., “michigan” in the first example of Table 1). (2)
“answer tokens” which contain the answer to the query (e.g., “motown” in the
second example of Table 1). The presence of such tokens may suggest the model
already “knows” the answer to the given question, either from pretraining or
from similar questions seen during training Lewis et al. (2021).

Figure 4:  The percentage of questions from the (entire) development set of NQ
with at least one query expansion token (i.e., a token that appears in the
question’s gold passage but not in the question itself) in the top-$k$ of the
question vocabulary projection $Q$, as a function of $k$. Stop words and
punctuation marks do not count as query expansion tokens.

Given these findings, we conjecture that the model “uses” these query
expansion tokens to introduce a semantic signal to the retrieval process.

Figure 5:  An analysis of token amnesia. We consider questions for which BM25
retrieves a correct passage (i.e., a passage that contains the answer) in its
top-5, and analyze what ranks were assigned to tokens shared by the question
and the passage in the passage vocabulary projection $P$ (left) and question
vocabulary projection $Q$ (right). We plot the maximal token rank as a
function of the rank assigned to the correct passage by DPR.

##  6 Token Amnesia

The analysis in Section 5 shows that vocabulary projections of passages (i.e.,
$P$) predict which of the input tokens are likely to appear in relevant
questions. However, in some cases these predictions utterly fail. For example,
in Figure 2 the token “michael” is missing from the top-$k$ of the passage
projection $P$. We refer to such cases as token amnesia. Here we ask, do these
failure in query prediction hurt retrieval?

Next, we demonstrate that token amnesia indeed correlates with well-known
failures of dense retrievers (§6.1). To overcome this issue, we suggest a
lexical enrichment procedure for dense representations (§6.2) and demonstrate
its effectiveness on downstream retrieval performance (§6.3).

###  6.1 Token Amnesia is Correlated with Retriever Failures

Dense retrievers have shown difficulties in out-of-domain settings Sciavolino
et al. (2021); Thakur et al. (2021), where even sparse models like BM25
significantly outperform them. We now offer an intuitive explanation to these
failures via token amnesia. We focus on setups where BM25 outperforms dense
models and ask: why do dense retrievers fail to model lexical overlap signals?
To answer this question, we consider subsets of NQ and EntityQs where BM25 is
able to retrieve the correct passage in its top-5 results. We focus on these
subsets as they contain significant lexical overlap between questions and
passages (by definition, as BM25 successfully retrieved the correct passage).
Let $q$ be a question and $p$ the passage retrieved by BM25 for $q$, and $Q$
and $P$ be their corresponding vocabulary projections for some dense
retriever. Also, let $\mathcal{T}\subseteq\mathcal{V}$ be the set of tokens
that appear in both $q$ and $p$ (excluding stop words). Figure 5 shows the
maximum (i.e., lowest) rank of tokens from $\mathcal{T}$ in the distributions
$P$ (left) and $Q$ (right) as a function of whether DPR is able to retrieve
this passage (i.e., the rank of $p$ in the retrieval results of DPR). Indeed,
the median max-rank over questions for which DPR succeeds to fetch $p$ in its
top-5 results (blue box) is much lower than that of questions for which DPR
fails to retrieve the passage (red box). As expected (due to the fact that
questions contain less tokens than passages), the ranks of shared tokens in
question projections $Q$ are much higher. However, the trend is present in $Q$
as well. Additional figures (for EntityQs; as well as median ranks instead of
max ranks) are given in App. C.

Overall, these findings indicate a correlation between token amnesia and
failures of DPR. Next, we introduce a method to address token amnesia in dense
retrievers, via lexical enrichment of dense representations.

###  6.2 Method: Lexical Enrichment

As suggested by the analysis in §6.1, dense retrievers have the tendency to
ignore some of their input tokens. We now leverage this insight to improve
these models. We refer to our method as lexical enrichment (LE) because it
enriches text encodings with specific lexical items.

Intuitively, a natural remedy to the “token amnesia” problem is to change the
retriever encoding such that it does include these tokens. For example, assume
the query $q$ is “Where was Michael Jack born?” and the corresponding passage
$p$ contains the text “Michael Jack was born in Folkestone, England”.
According to Figure 2, the token “michael” is ranked relatively low in $P$,
and DPR fails to retrieve the correct passage $p$. We would like to modify the
passage representation $\bm{e}_{p}$ and get an enriched version
$\bm{e}^{\prime}_{p}$ that does have this token in its top-$k$ projected
tokens, while keeping most of the other projected tokens intact. This is our
goal in LE, and we next describe the approach. We focus on enrichment of
passage representations, as query enrichment works similarly. We first explain
how to enrich representations with a single token, and then extend the process
to multiple tokens.

#### Single-Token Enrichment

Assume we want to enrich a passage representation $\bm{e}_{p}$ with a token
$t$ (e.g., $t=``$michael” in the above example). If there were no other words
in the passage, we’d simply want to find an embedding such that feeding it
into the MLM would produce $t$ as the top token.333Note that feeding the token
input embedding $\bm{v}_{t}$ does not necessarily produce $t$ as the top
token, as the MLM head applies a non-linear function $g$ (Eq. 1). We refer to
this embedding as the single-token enrichment of $t$, denote it by
$\bm{s}_{t}$ and define it as:444This is equivalent to the cross-entropy loss
between a one-hot vector on $t$ and the output distribution
$MLM(\hat{\bm{s}})$.

$\bm{s}_{t}=\operatorname*{arg\,max}_{\hat{\bm{s}}}\ \ \log\text{MLM-
Head}(\hat{\bm{s}})[t]$ (4)

In order to approximately solve the optimization problem in Eq. 4 for each $t$
in the vocabulary, we use Adam with a learning rate of 0.01.555For S-MPNet, we
used a learning rate of $10^{-3}$. We stop when a (cross-entropy) loss
threshold of 0.1 is reached for all tokens. We then apply whitening Jung et
al. (2022), which was proven effective for dense retrieval.

Model | $\lambda$ |  BEIR |  MTEB | EntityQs | TriviaQA |  WQ |  TREC | SQuAD
---|---|---|---|---|---|---|---|---
nDCG@10 | Top-20 retrieval accuracy
BM25 | - | 42.9 | 42.3 | 71.4 | 76.4 | 62.4 | 81.1 | 71.2
BM25 (BERT/MPNet Tokens) | - | 41.6 | 41.7 | 66.2 | 75.8 | 62.1 | 79.3 | 70.0
DPR | - | 21.4 | 22.4 | 49.7 | 69.0 | 68.8 | 85.9 | 48.9
DPR + LE | 5.0 | 26.4 | 27.6 | 65.4 | 75.3 | 73.2 | 87.9 | 59.7
S-MPNet | - | 43.1 | 44.6 | 57.6 | 77.6 | 73.9 | 90.2 | 65.5
S-MPNet + LE | 0.5 | 44.1 | 45.7 | 68.5 | 78.9 | 74.5 | 90.4 | 69.0
Spider | - | 27.4 | 26.4 | 66.3 | 75.8 | 65.9 | 82.6 | 61.0
Spider + LE | 3.0 | 29.5 | 28.8 | 68.9 | 76.3 | 70.2 | 83.4 | 62.8
Table 3: Retrieval results on BEIR, the retrieval cluster of MTEB and five
open-domain QA datasets. LE stands for lexical enrichment (our method; §6.2),
that enriches query and passage representation with lexical information.
$\lambda$ is defined in Eq. 5. BM25 (BERT Vocabulary) refers to a model that
operates over tokens from BERT’s vocabulary, rather than words. For each model
and dataset, we compare the enriched (LE) model with the original, and mark in
bold the better one from the two. We underline the best overall model for each
dataset. Results for each of the BEIR datasets are given in Table 9.
Top-$\\{1,5,100\\}$ accuracy results are given in Tables 6, 7 & 8.

#### Multi-Token Enrichment

Now suppose we have an input $x$ (either a question or a passage) and we’d
like to enrich its representation with its tokens $x=[x_{1},..,x_{n}]$, such
that rare tokens are given higher weights than frequent ones (as in BM25).
Then, we simply take its original representation $\bm{e}_{x}$ and add to it a
weighted sum of the single-token enrichments (Eq. 4). Namely, we define:

$\begin{split}\bm{e}_{x}^{\text{lex}}&=\frac{1}{n}\sum_{i=1}^{n}w_{x_{i}}\bm{s}_{x_{i}}\\\
\bm{e}^{\prime}_{x}&=\bm{e}_{x}+\lambda\cdot\frac{\bm{e}_{x}^{\text{lex}}}{||\bm{e}_{x}^{\text{lex}}||}\end{split}$
(5)

Here $\lambda$ is a hyper-parameter chosen via cross validation. We use the
inverse document frequency Sparck Jones (1972) of tokens as their weights:
$w_{x_{i}}=\text{IDF}(x_{i})$. The relevance score is then defined on the
enriched representations.

###  6.3 Results

Our experiments demonstrate the effectiveness of our method for multiple
models, especially in zero-shot settings. Table 3 shows the results of several
models with and without our enrichment method, LE. Additional results are
given in App. D. The results demonstrate the effectiveness of LE when added to
all baseline models. Importantly, our method improves the performance of
S-MPNet—the best base-sized model on the MTEB benchmark to date Muennighoff et
al. (2022)—on MTEB and BEIR by 1.1% and 1.0%, respectively. When considering
EntityQs (on which dense retrievers are known to struggle), we observe
significant gains across all models, and S-MPNet and Spider obtain higher
accuracy than BM25 that operates on the same textual units (i.e., BM25 with
BERT vocabulary). This finding indicates that they are able to integrate
semantic information (from the original representation) with lexical signals.
Yet, vanilla BM25 is still better than LE models on EntityQs and SQuAD, which
prompts further work on how to incorporate lexical signals in dense retrieval.
Overall, it is evident that LE improves retrieval accuracy compared to
baseline models for all models and datasets (i.e., zero-shot setting).

Method |  | NQ (Dev Set) |  | EntityQs (Dev Set)
---|---|---|---|---
| Top-1 | Top-5 | Top-20 | Top-100 |  | Top-1 | Top-5 | Top-20 | Top-100
DPR |  | 44.9 | 66.8 | 78.1 | 85.0 |  | 24.0 | 38.4 | 50.4 | 63.5
DPR + LE |  | 44.4 | 67.5 | 79.4 | 86.0 |  | 38.3 | 54.0 | 65.2 | 76.1
No IDF |  | 45.1 | 67.3 | 78.5 | 85.4 |  | 32.0 | 46.4 | 57.7 | 69.6
BERT embedding matrix |  | 44.8 | 67.6 | 79.1 | 85.6 |  | 34.6 | 50.3 | 61.8 | 72.8
No whitening |  | 44.1 | 66.3 | 78.7 | 85.2 |  | 34.6 | 49.7 | 61.4 | 72.9
No $\ell_{2}$ normalization |  | 43.9 | 66.8 | 79.2 | 86.0 |  | 35.5 | 51.3 | 63.0 | 74.6
Table 4: Ablation study on the development set of Natrual Questions and Entity
Questions. DPR + LE is our lexical enrichment method applied on DPR. No IDF
removes the IDF weights in Eq. 5 (i.e., mean pooling). BERT embedding matrix
replaces single-token enrichment $\bm{s}_{t}$ as defined in Eq. 4 with the
static token embeddings of BERT, $\bm{v}_{t}$ (Eq. 1). No whitening removes
whitening transformation. No $\ell_{2}$ normalization removes the
normalization of $\bm{e}_{x}^{\text{lex}}$.

###  6.4 Ablation Study

We carry an ablation study to test our design choices from §6.2. We evaluate
four elements of our method: (1) The use of IDF to highlight rare tokens, (2)
Our approach for deriving single-token representations, (3) The use of
whitening, and (4) The use of unit normalization.

#### IDF

In our method, we create lexical representations of questions and passages,
$\bm{e}_{x}^{\text{lex}}$. These lexical representations are the average of
token embeddings, each multiplied by its token’s IDF. We validate that IDF is
indeed necessary – Table 4 demonstrates that setting $w_{x_{i}}=1$ in Eq. 5
leads to a significant degradation in performance on EntityQs. For example,
top-20 retrieval accuracy drops from 65.2% to 57.7%.

#### Single-Token Enrichment

Eq. 4 defines our single-token enrichment: for each item in the vocabulary
$v\in\mathcal{V}$, we find an embedding which gives a one-hot vector peaked at
$v$ when fed to the MLM head. We confirm that this is necessary by replacing
Eq. 4 with the static embeddings of the pretrained model (e.g., BERT in the
case of DPR). We find that our approach significantly improves over BERT’s
embeddings on EntityQs (e.g., the margin in top-20 accuracy is 3.4%).

#### Whitening & Normalization

Last, we experiment with removing the whitening and $\ell_{2}$ normalization.
It is evident that they are both necessary, as removing either of them causes
a dramatic drop in performance (3.8% and 2.2% in top-20 accuracy on EntityQs,
respectively).

##  7 Related Work

Projecting representations and model parameters to the vocabulary space has
been studied previously mainly in the context of language models. The approach
was initially explored by nostalgebraist (2020). Geva et al. (2021) showed
that feed-forward layers in transformers can be regarded as key-value
memories, where the value vectors induce distributions over the vocabulary.
Geva et al. (2022) view the token representations themselves as inducing such
distributions, with feed-forward layers “updating” them. Dar et al. (2022)
suggest to project all transformer parameters to the vocabulary space. Dense
retrieval models, however, do not have any language modeling objective during
fine-tuning, yet we show that their representations can still be projected to
the vocabulary.

Despite the wide success of dense retrievers recently, interpreting their
representations remains under-explored. MacAvaney et al. (2022) analyze neural
retrieval models (not only dense retrievers) via diagnostic probes, testing
characteristics like sensitivity to paraphrases, styles and factuality.
Adolphs et al. (2022) decode the query representations of neural retrievers
using a T5 decoder, and show how to “move” in representation space to decode
better queries for retrieval.

Language models (and specifically MLMs) have been used for sparse retrieval in
the context of term-weighting and lexical expansion. For example, Bai et al.
(2020) and Formal et al. (2021) learn such functions over BERT’s vocabulary
space. We differ by showing that dense retrievers implicitly operate in that
space as well. Thus, these approaches may prove effective for dense models as
well. While we focus in this work on dense retrievers based on encoder-only
models, our framework is easily extendable for retrievers based on
autoregressive decoder-only (i.e., left-to-right) models like GPT Radford et
al. (2019); Brown et al. (2020), e.g., Neelakantan et al. (2022) and
Muennighoff (2022).

##  8 Conclusion

In this work, we explore projecting query and passage representations obtained
by dense retrieval to the vocabulary space. We show that these projections
facilitate a better understanding of the mechanisms underlying dense
retrieval, as well as their failures. We also demonstrate how projections can
help improve these models. This understanding is likely to help in improving
retrievers, as our lexical enrichment approach demonstrates.

## Limitations

We point to several limitations of our work. First, our work considers a
popular family of models referred to as “dense retrievers”, but other
approaches for retrieval include sparse retrievers Robertson and Zaragoza
(2009); Bai et al. (2020); Formal et al. (2021), generative retrievers Tay et
al. (2022); Bevilacqua et al. (2022), late-interaction models Khattab and
Zaharia (2020), inter alia. While our work draws interesting connections
between dense and sparse retrieval, our main focus is on understanding and
improving dense models. Second, all three dense models we analyze are
bidirectional and were trained in a contrastive fashion. While most dense
retrievers indeed satisfy these properties, there are works that suggested
other approaches, both in terms of other architectures Muennighoff (2022);
Neelakantan et al. (2022); Ni et al. (2022) and other training frameworks
Lewis et al. (2020); Izacard et al. (2022b). Last, while our work introduces
new ways to interpret and analyze dense retrieval models, we believe our work
is the tip of the iceberg, and there is still much work to be done in order to
gain a full understanding of these models.

## Ethics Statement

Retrieval systems have the potential to mitigate serious problems caused by
language models, like factual inaccuracies. However, retrieval failures may
lead to undesirable behavior of downstream models, like wrong answers in QA or
incorrect generations for other tasks. Also, since retrieval models are based
on pretrained language models, they may suffer from similar biases.

## Acknowledgements

We thank Ori Yoran, Yoav Levine, Yuval Kirstain, Mor Geva and the anonymous
reviewers for their valuable feedback. This project was funded by the European
Research Council (ERC) under the European Unions Horizon 2020 research and
innovation programme (grant ERC HOLI 819080), the Blavatnik Fund, the Alon
Scholarship, the Yandex Initiative for Machine Learning, Intel Corporation,
ISRAEL SCIENCE FOUNDATION (grant No. 448/20), Open Philanthropy, and an
Azrieli Foundation Early Career Faculty Fellowship.

## References

  * Adolphs et al. (2022) Leonard Adolphs, Michelle Chen Huebscher, Christian Buck, Sertan Girgin, Olivier Bachem, Massimiliano Ciaramita, and Thomas Hofmann. 2022.  Decoding a neural retriever’s latent space for query suggestion.
  * Bai et al. (2020) Yang Bai, Xiaoguang Li, Gang Wang, Chaoliang Zhang, Lifeng Shang, Jun Xu, Zhaowei Wang, Fangshan Wang, and Qun Liu. 2020.  SparTerm: Learning term-based sparse representation for fast text retrieval.
  * Baudiš and Šedivý (2015) Petr Baudiš and Jan Šedivý. 2015.  Modeling of the question answering task in the YodaQA system.  In _Proceedings of the 6th International Conference on Experimental IR Meets Multilinguality, Multimodality, and Interaction - Volume 9283_ , CLEF’15, page 222–228, Berlin, Heidelberg. Springer-Verlag.
  * Berant et al. (2013) Jonathan Berant, Andrew Chou, Roy Frostig, and Percy Liang. 2013.  Semantic parsing on Freebase from question-answer pairs.  In _Proceedings of the 2013 Conference on Empirical Methods in Natural Language Processing_ , pages 1533–1544, Seattle, Washington, USA. Association for Computational Linguistics.
  * Bevilacqua et al. (2022) Michele Bevilacqua, Giuseppe Ottaviano, Patrick Lewis, Wen-tau Yih, Sebastian Riedel, and Fabio Petroni. 2022.  Autoregressive search engines: Generating substrings as document identifiers.
  * Brown et al. (2020) Tom B. Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, Sandhini Agarwal, Ariel Herbert-Voss, Gretchen Krueger, Tom Henighan, Rewon Child, Aditya Ramesh, Daniel M. Ziegler, Jeffrey Wu, Clemens Winter, Christopher Hesse, Mark Chen, Eric Sigler, Mateusz Litwin, Scott Gray, Benjamin Chess, Jack Clark, Christopher Berner, Sam McCandlish, Alec Radford, Ilya Sutskever, and Dario Amodei. 2020.  Language models are few-shot learners.  In _Advances in Neural Information Processing Systems_.
  * Chang et al. (2020) Wei-Cheng Chang, Felix X. Yu, Yin-Wen Chang, Yiming Yang, and Sanjiv Kumar. 2020\.  Pre-training tasks for embedding-based large-scale retrieval.  In _8th International Conference on Learning Representations, ICLR 2020_.
  * Chen et al. (2022) Xilun Chen, Kushal Lakhotia, Barlas Oğuz, Anchit Gupta, Patrick Lewis, Stan Peshterliev, Yashar Mehdad, Sonal Gupta, and Wen-tau Yih. 2022.  Salient phrase aware dense retrieval: Can a dense retriever imitate a sparse one?
  * Dar et al. (2022) Guy Dar, Mor Geva, Ankit Gupta, and Jonathan Berant. 2022.  Analyzing transformers in embedding space.
  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019.  BERT: Pre-training of deep bidirectional transformers for language understanding.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , pages 4171–4186, Minneapolis, Minnesota. Association for Computational Linguistics.
  * Formal et al. (2021) Thibault Formal, Benjamin Piwowarski, and Stéphane Clinchant. 2021.  SPLADE: Sparse lexical and expansion model for first stage ranking.  In _Proceedings of the 44th International ACM SIGIR Conference on Research and Development in Information Retrieval_ , SIGIR ’21, page 2288–2292, New York, NY, USA. Association for Computing Machinery.
  * Geva et al. (2022) Mor Geva, Avi Caciularu, Kevin Ro Wang, and Yoav Goldberg. 2022.  Transformer feed-forward layers build predictions by promoting concepts in the vocabulary space.
  * Geva et al. (2021) Mor Geva, Roei Schuster, Jonathan Berant, and Omer Levy. 2021.  Transformer feed-forward layers are key-value memories.  In _Proceedings of the 2021 Conference on Empirical Methods in Natural Language Processing_ , pages 5484–5495, Online and Punta Cana, Dominican Republic. Association for Computational Linguistics.
  * Izacard et al. (2022a) Gautier Izacard, Mathilde Caron, Lucas Hosseini, Sebastian Riedel, Piotr Bojanowski, Armand Joulin, and Edouard Grave. 2022a.  Unsupervised dense information retrieval with contrastive learning.  _Transactions on Machine Learning Research_.
  * Izacard et al. (2022b) Gautier Izacard, Patrick Lewis, Maria Lomeli, Lucas Hosseini, Fabio Petroni, Timo Schick, Jane Dwivedi-Yu, Armand Joulin, Sebastian Riedel, and Edouard Grave. 2022b.  Atlas: Few-shot learning with retrieval augmented language models.
  * Johnson et al. (2021) Jeff Johnson, Matthijs Douze, and Hervé Jégou. 2021.  Billion-scale similarity search with GPUs.  _IEEE Transactions on Big Data_ , 7(3):535–547.
  * Joshi et al. (2017) Mandar Joshi, Eunsol Choi, Daniel Weld, and Luke Zettlemoyer. 2017.  TriviaQA: A large scale distantly supervised challenge dataset for reading comprehension.  In _Proceedings of the 55th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 1601–1611, Vancouver, Canada. Association for Computational Linguistics.
  * Jung et al. (2022) Euna Jung, Jungwon Park, Jaekeol Choi, Sungyoon Kim, and Wonjong Rhee. 2022.  Isotropic representation can improve dense retrieval.
  * Karpukhin et al. (2020) Vladimir Karpukhin, Barlas Oguz, Sewon Min, Patrick Lewis, Ledell Wu, Sergey Edunov, Danqi Chen, and Wen-tau Yih. 2020.  Dense passage retrieval for open-domain question answering.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pages 6769–6781, Online. Association for Computational Linguistics.
  * Khattab and Zaharia (2020) Omar Khattab and Matei Zaharia. 2020.  ColBERT: Efficient and effective passage search via contextualized late interaction over BERT.  In _Proceedings of the 43rd International ACM SIGIR Conference on Research and Development in Information Retrieval_ , SIGIR ’20, page 39–48, New York, NY, USA. Association for Computing Machinery.
  * Kwiatkowski et al. (2019) Tom Kwiatkowski, Jennimaria Palomaki, Olivia Redfield, Michael Collins, Ankur Parikh, Chris Alberti, Danielle Epstein, Illia Polosukhin, Jacob Devlin, Kenton Lee, Kristina Toutanova, Llion Jones, Matthew Kelcey, Ming-Wei Chang, Andrew M. Dai, Jakob Uszkoreit, Quoc Le, and Slav Petrov. 2019.  Natural questions: A benchmark for question answering research.  _Transactions of the Association for Computational Linguistics_ , 7:452–466.
  * Lee et al. (2019) Kenton Lee, Ming-Wei Chang, and Kristina Toutanova. 2019.  Latent retrieval for weakly supervised open domain question answering.  In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_ , pages 6086–6096, Florence, Italy. Association for Computational Linguistics.
  * Lewis et al. (2020) Patrick Lewis, Ethan Perez, Aleksandra Piktus, Fabio Petroni, Vladimir Karpukhin, Naman Goyal, Heinrich Küttler, Mike Lewis, Wen-tau Yih, Tim Rocktäschel, Sebastian Riedel, and Douwe Kiela. 2020.  Retrieval-augmented generation for knowledge-intensive nlp tasks.  In _Advances in Neural Information Processing Systems_ , volume 33, pages 9459–9474. Curran Associates, Inc.
  * Lewis et al. (2021) Patrick Lewis, Pontus Stenetorp, and Sebastian Riedel. 2021.  Question and answer test-train overlap in open-domain question answering datasets.  In _Proceedings of the 16th Conference of the European Chapter of the Association for Computational Linguistics: Main Volume_ , pages 1000–1008, Online. Association for Computational Linguistics.
  * Lin et al. (2021) Jimmy Lin, Xueguang Ma, Sheng-Chieh Lin, Jheng-Hong Yang, Ronak Pradeep, and Rodrigo Nogueira. 2021.  Pyserini: A python toolkit for reproducible information retrieval research with sparse and dense representations.  In _Proceedings of the 44th International ACM SIGIR Conference on Research and Development in Information Retrieval_ , SIGIR ’21, page 2356–2362, New York, NY, USA. Association for Computing Machinery.
  * Liu et al. (2019) Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  RoBERTa: A robustly optimized bert pretraining approach.
  * MacAvaney et al. (2022) Sean MacAvaney, Sergey Feldman, Nazli Goharian, Doug Downey, and Arman Cohan. 2022\.  ABNIRML: Analyzing the behavior of neural IR models.  _Transactions of the Association for Computational Linguistics_ , 10:224–239.
  * Mao et al. (2021) Yuning Mao, Pengcheng He, Xiaodong Liu, Yelong Shen, Jianfeng Gao, Jiawei Han, and Weizhu Chen. 2021.  Generation-augmented retrieval for open-domain question answering.  In _Proceedings of the 59th Annual Meeting of the Association for Computational Linguistics and the 11th International Joint Conference on Natural Language Processing (Volume 1: Long Papers)_ , pages 4089–4100, Online. Association for Computational Linguistics.
  * Muennighoff (2022) Niklas Muennighoff. 2022.  SGPT: GPT sentence embeddings for semantic search.
  * Muennighoff et al. (2022) Niklas Muennighoff, Nouamane Tazi, Loïc Magne, and Nils Reimers. 2022.  MTEB: Massive text embedding benchmark.
  * Neelakantan et al. (2022) Arvind Neelakantan, Tao Xu, Raul Puri, Alec Radford, Jesse Michael Han, Jerry Tworek, Qiming Yuan, Nikolas Tezak, Jong Wook Kim, Chris Hallacy, Johannes Heidecke, Pranav Shyam, Boris Power, Tyna Eloundou Nekoul, Girish Sastry, Gretchen Krueger, David Schnurr, Felipe Petroski Such, Kenny Hsu, Madeleine Thompson, Tabarak Khan, Toki Sherbakov, Joanne Jang, Peter Welinder, and Lilian Weng. 2022.  Text and code embeddings by contrastive pre-training.
  * Ni et al. (2022) Jianmo Ni, Chen Qu, Jing Lu, Zhuyun Dai, Gustavo Hernández Ábrego, Ji Ma, Vincent Y. Zhao, Yi Luan, Keith B. Hall, Ming-Wei Chang, and Yinfei Yang. 2022\.  Large dual encoders are generalizable retrievers.
  * nostalgebraist (2020) nostalgebraist. 2020.  interpreting gpt: the logit lens.
  * Qu et al. (2021) Yingqi Qu, Yuchen Ding, Jing Liu, Kai Liu, Ruiyang Ren, Wayne Xin Zhao, Daxiang Dong, Hua Wu, and Haifeng Wang. 2021.  RocketQA: An optimized training approach to dense passage retrieval for open-domain question answering.  In _Proceedings of the 2021 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_ , pages 5835–5847, Online. Association for Computational Linguistics.
  * Radford et al. (2019) Alec Radford, Jeff Wu, Rewon Child, David Luan, Dario Amodei, and Ilya Sutskever. 2019.  Language models are unsupervised multitask learners.
  * Rajpurkar et al. (2016) Pranav Rajpurkar, Jian Zhang, Konstantin Lopyrev, and Percy Liang. 2016.  SQuAD: 100,000+ questions for machine comprehension of text.  In _Proceedings of the 2016 Conference on Empirical Methods in Natural Language Processing_ , pages 2383–2392, Austin, Texas. Association for Computational Linguistics.
  * Ram et al. (2022) Ori Ram, Gal Shachaf, Omer Levy, Jonathan Berant, and Amir Globerson. 2022.  Learning to retrieve passages without supervision.  In _Proceedings of the 2022 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_ , pages 2687–2700, Seattle, United States. Association for Computational Linguistics.
  * Reddy et al. (2021) Revanth Gangi Reddy, Vikas Yadav, Md Arafat Sultan, Martin Franz, Vittorio Castelli, Heng Ji, and Avirup Sil. 2021.  Towards robust neural retrieval models with synthetic pre-training.
  * Reimers and Gurevych (2019) Nils Reimers and Iryna Gurevych. 2019.  Sentence-BERT: Sentence embeddings using Siamese BERT-networks.  In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_ , pages 3982–3992, Hong Kong, China. Association for Computational Linguistics.
  * Robertson and Zaragoza (2009) Stephen Robertson and Hugo Zaragoza. 2009.  The probabilistic relevance framework: BM25 and beyond.  _Found. Trends Inf. Retr._ , 3(4):333–389.
  * Rocchio (1971) Joseph Rocchio. 1971.  Relevance feedback in information retrieval.  _The SMART retrieval system: experiments in automatic document processing_ , pages 313–323.
  * Sciavolino et al. (2021) Christopher Sciavolino, Zexuan Zhong, Jinhyuk Lee, and Danqi Chen. 2021.  Simple entity-centric questions challenge dense retrievers.  In _Proceedings of the 2021 Conference on Empirical Methods in Natural Language Processing_ , pages 6138–6148, Online and Punta Cana, Dominican Republic. Association for Computational Linguistics.
  * Song et al. (2020) Kaitao Song, Xu Tan, Tao Qin, Jianfeng Lu, and Tie-Yan Liu. 2020.  MPNet: Masked and permuted pre-training for language understanding.  In _Advances in Neural Information Processing Systems_ , volume 33, pages 16857–16867. Curran Associates, Inc.
  * Sparck Jones (1972) Karen Sparck Jones. 1972.  A statistical interpretation of term specificity and its application in retrieval.  In _Journal of Documentation_ , volume 28 no. 1, pages 11–21.
  * Tay et al. (2022) Yi Tay, Vinh Q. Tran, Mostafa Dehghani, Jianmo Ni, Dara Bahri, Harsh Mehta, Zhen Qin, Kai Hui, Zhe Zhao, Jai Gupta, Tal Schuster, William W. Cohen, and Donald Metzler. 2022.  Transformer memory as a differentiable search index.  In _Advances in Neural Information Processing Systems_.
  * Thakur et al. (2021) Nandan Thakur, Nils Reimers, Andreas Rücklé, Abhishek Srivastava, and Iryna Gurevych. 2021.  BEIR: A heterogeneous benchmark for zero-shot evaluation of information retrieval models.  In _Proceedings of the Neural Information Processing Systems Track on Datasets and Benchmarks_ , volume 1.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Ł ukasz Kaiser, and Illia Polosukhin. 2017.  Attention is all you need.  In _Advances in Neural Information Processing Systems 30_ , pages 5998–6008.
  * Voorhees (1994) Ellen M. Voorhees. 1994.  Query expansion using lexical-semantic relations.  In _Proceedings of the 17th Annual International ACM SIGIR Conference on Research and Development in Information Retrieval_ , SIGIR ’94, page 61–69, Berlin, Heidelberg. Springer-Verlag.
  * Wolf et al. (2020) Thomas Wolf, Lysandre Debut, Victor Sanh, Julien Chaumond, Clement Delangue, Anthony Moi, Pierric Cistac, Tim Rault, Remi Louf, Morgan Funtowicz, Joe Davison, Sam Shleifer, Patrick von Platen, Clara Ma, Yacine Jernite, Julien Plu, Canwen Xu, Teven Le Scao, Sylvain Gugger, Mariama Drame, Quentin Lhoest, and Alexander Rush. 2020.  Transformers: State-of-the-art natural language processing.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing: System Demonstrations_ , pages 38–45, Online. Association for Computational Linguistics.
  * Zhao and Callan (2012) Le Zhao and Jamie Callan. 2012.  Automatic term mismatch diagnosis for selective query expansion.  In _Proceedings of the 35th International ACM SIGIR Conference on Research and Development in Information Retrieval_ , SIGIR ’12, page 515–524, New York, NY, USA. Association for Computing Machinery.

##  Appendix A Models: Further Details

#### DPR

Karpukhin et al. (2020) is a dense retriever that was trained on Natural
Questions Kwiatkowski et al. (2019). It was initialized from BERT-base Devlin
et al. (2019). Thus, we use the public pretrained MLM head of BERT-base to
project DPR representations.

#### BERT

Devlin et al. (2019) We use BERT for dense retrieval, mainly as a baseline for
DPR, as DPR was initialized from BERT. This allows us to track where behaviors
we observe stem from: pretraining or retrieval fine-tuning. We use both CLS
and mean pooling for BERT.

#### S-MPNet

is a supervised model trained for Sentence Transformers Reimers and Gurevych
(2019) using many available datasets for retrieval, sentence similarity, inter
alia. It uses cosine similarity, rather than dot product, for relevance
scores. It was initialized from MPNet-base Song et al. (2020), and thus we use
this model’s MLM head.

#### Spider

Ram et al. (2022) is an unsupervised dense retriever trained using the
recurring span retrieval pretraining task. It was also initialized from BERT-
base, and we therefore use the same MLM head for projection as the one used
for DPR.

#### BM25

Robertson and Zaragoza (2009) is a lexical model based on tf-idf. We use two
variants of BM25: (1) vanilla BM25, and (2) BM25 over BERT/MPNet tokens (e.g.,
“Reba” $\rightarrow$ “re ##ba”).666BERT and MPNet use essentially the same
vocabulary, up to special tokens. We consider this option to understand
whether the advantages of BM25 stem from its use of different word units from
the transformer models.

##  Appendix B Analysis: Further Results

Figure 6 gives an analysis of the top-$k$ tokens in the question projection
$Q$ and passage projection $P$.

##  Appendix C Token Amnesia: Further results

Figure 7 gives further analyses of token amnesia: It contains the results for
EntityQuestions, as well as analysis of median ranks in addition to max ranks
(complements Figure 5).

##  Appendix D Lexical Enrichment: Further Results

(a) Question projection $Q$ (b) Passage projection $P$

Figure 6: An analysis of the top-$k$ tokens in the vocabulary projection $Q$ (a) for questions from the development set of NQ and $P$ (b) for their corresponding gold passage of DPR. Specifically, we analyze what percentage of these top-$k$ tokens are present in the question and/or the passage for $k\in\\{1,5,20\\}$.  Dataset | License | Test Ex.
---|---|---
Natural Questions | Apache-2.0 | 3,610
TriviaQA | Apache-2.0 | 11,313
WebQuestions | CC BY 4.0 | 2,032
CuratedTREC | - | 694
SQuAD | CC BY-SA 4.0 | 10,570
EntityQs | MIT | 22,075
Table 5: The license and number of test example in each of the datasets used
in the paper.

Table 9 gives the results of our method on the BEIR and MTEB benchmarks for
all 19 datasets (complements Table 3). Table 6, Table 7 and Table 8 give the
zero-shot results for $k\in\\{1,5,100\\}$, respectively (complement Table 3).

##  Appendix E Dataset Statistics & Licenses

Table 5 details the license and number of test example for each of the six
open-domain datasets used in our work. For the BEIR benchmark, we refer the
reader to Thakur et al. (2021) for number of examples and license of each of
their datasets.

##  Appendix F Computational Resources

Our method (LE) does not involve training models at all. Our computational
resources have been used to evaluate LE on the BEIR benchmark, i.e., computing
passage embeddings for each corpus and each model. We used eight Quadro RTX
8000 GPUs. Each experiment took several hours.

(a) Max rank among shared tokens, EntityQuestions (b) Median rank among shared
tokens, Natural Questions (c) Median rank among shared tokens, EntityQuestions

Figure 7:  Further analysis of token amnesia (complementary to Figure 5). We consider questions for which BM25 retrieves a correct passage (i.e., a passage that contains the answer) in its top-5, and analyze what ranks were assigned to tokens shared by the question and the passage in the passage vocabulary projection $P$ (left) and question vocabulary projection $Q$ (right). We plot the max and median token rank as a function of the rank assigned to the correct passage by DPR, for Natural Questions (NQ) and EntityQuestions (EQ).  Model | EntityQs | TriviaQA |  WQ |  TREC | SQuAD
---|---|---|---|---|---
BM25 | 43.5 | 46.3 | 18.9 | 34.6 | 36.7
BM25 (BERT/MPNet Vocabulary) | 37.6 | 45.4 | 19.2 | 33.0 | 35.6
DPR | 24.3 | 37.3 | 30.5 | 51.3 | 16.0
DPR + LE | 38.3 | 45.8 | 35.0 | 54.6 | 22.8
S-MPNet | 22.7 | 42.9 | 30.9 | 51.0 | 25.8
S-MPNet + LE | 37.3 | 47.3 | 37.1 | 54.0 | 30.0
Spider | 35.0 | 41.7 | 22.3 | 38.2 | 22.2
Spider + LE | 40.7 | 43.7 | 27.8 | 43.2 | 23.5
Table 6: Top-1 retrieval accuracy in a “zero-shot” setting (i.e., datasets were not used for model training), complementary to Table 3. LE stands for lexical enrichment (our method; §6.2), that enriches query and passage representation with lexical information. BM25 (BERT Vocabulary) refers to a model that operates over tokens from BERT’s vocabulary, rather than words. For each model and dataset, we compare the enriched (LE) model with the original, and mark in bold the better one from the two. We underline the best overall model for each dataset.  Model | EntityQs | TriviaQA |  WQ |  TREC | SQuAD
---|---|---|---|---|---
BM25 | 61.0 | 66.3 | 41.8 | 64.6 | 57.5
BM25 (BERT/MPNet Vocabulary) | 55.1 | 65.6 | 42.3 | 62.5 | 56.1
DPR | 38.1 | 57.0 | 52.7 | 74.1 | 33.4
DPR + LE | 53.8 | 64.8 | 57.7 | 79.5 | 42.3
S-MPNet | 42.7 | 66.1 | 58.8 | 79.7 | 49.5
S-MPNet + LE | 56.8 | 68.5 | 61.6 | 81.4 | 53.2
Spider | 54.5 | 63.6 | 46.8 | 65.9 | 43.6
Spider + LE | 58.0 | 64.4 | 52.2 | 70.0 | 44.9
Table 7: Top-5 retrieval accuracy in a “zero-shot” setting (i.e., datasets were not used for model training), complementary to Table 3. LE stands for lexical enrichment (our method; §6.2), that enriches query and passage representation with lexical information. BM25 (BERT Vocabulary) refers to a model that operates over tokens from BERT’s vocabulary, rather than words. For each model and dataset, we compare the enriched (LE) model with the original, and mark in bold the better one from the two. We underline the best overall model for each dataset.  Model | EntityQs | TriviaQA |  WQ |  TREC | SQuAD
---|---|---|---|---|---
BM25 | 80.0 | 83.2 | 75.5 | 90.3 | 82.0
BM25 (BERT/MPNet Vocabulary) | 76.6 | 83.0 | 76.0 | 90.5 | 81.1
DPR | 63.2 | 78.7 | 78.3 | 92.1 | 65.1
DPR + LE | 76.1 | 82.9 | 82.1 | 93.5 | 74.0
S-MPNet | 71.7 | 84.8 | 83.0 | 95.1 | 78.4
S-MPNet + LE | 78.6 | 85.1 | 83.8 | 95.0 | 80.7
Spider | 77.4 | 83.5 | 79.7 | 92.8 | 76.0
Spider + LE | 78.9 | 83.8 | 81.5 | 92.2 | 77.8
Table 8: Top-100 retrieval accuracy in a “zero-shot” setting (i.e., datasets were not used for model training), complementary to Table 3. LE stands for lexical enrichment (our method; §6.2), that enriches query and passage representation with lexical information. BM25 (BERT Vocabulary) refers to a model that operates over tokens from BERT’s vocabulary, rather than words. For each model and dataset, we compare the enriched (LE) model with the original, and mark in bold the better one from the two. We underline the best overall model for each dataset.  Dataset | DPR |  | Spider |  | S-MPNet
---|---|---|---|---|---
Original |  \+ LE |  | Original |  \+ LE |  | Original |  \+ LE
MS MARCO | 18.4 | 20.9 |  | 14.6 | 16.2 |  | 40.0 | 40.3
TREC-COVID | 22.2 | 30.8 |  | 30.5 | 32.0 |  | 51.0 | 51.3
NFCorpus | 15.7 | 19.0 |  | 27.4 | 26.2 |  | 33.4 | 33.6
NQ | 51.3 | 49.8 |  | 12.6 | 17.0 |  | 52.2 | 52.8
HotpotQA | 32.6 | 37.7 |  | 40.4 | 43.1 |  | 45.2 | 48.3
FiQA-2018 | 10.5 | 13.0 |  |  1.0 | 11.2 |  | 49.3 | 49.8
ArguAna | 10.8 | 14.1 |  | 31.2 | 31.0 |  | 39.6 | 49.2
Touché-2020 | 13.1 | 15.8 |  |  4.2 |  6.4 |  | 21.0 | 21.5
CQADupStack | 12.7 | 18.0 |  | 21.3 | 21.7 |  | 44.6 | 44.7
Quora | 16.8 | 42.4 |  | 73.0 | 75.6 |  | 87.0 | 87.3
DBPedia | 26.9 | 28.5 |  | 20.0 | 22.3 |  | 34.1 | 34.8
SCIDOCS |  7.4 | 10.1 |  | 13.1 | 12.8 |  | 23.6 | 23.5
FEVER | 52.7 | 54.7 |  | 30.2 | 34.3 |  | 59.0 | 60.0
Climate-FEVER | 18.2 | 22.9 |  | 12.4 | 22.4 |  | 23.1 | 23.6
SciFact | 26.9 | 36.1 |  | 63.6 | 59.8 |  | 65.2 | 65.3
BioASQ | 11.6 | 17.6 |  | 21.0 | 22.3 |  | 21.5 | 22.3
Signal-1M (RT) | 13.6 | 21.1 |  | 25.3 | 26.1 |  | 24.9 | 25.3
TREC-NEWS | 19.1 | 21.3 |  | 29.3 | 31.3 |  | 50.7 | 50.7
Robust04 | 22.4 | 22.7 |  | 36.4 | 35.9 |  | 50.0 | 50.0
Avg. (MTEB: Retrieval) | 22.4 | 27.6 |  | 26.4 | 28.8 |  | 44.6 | 45.7
Avg. (BEIR) | 21.4 | 26.4 |  | 27.4 | 29.5 |  | 43.1 | 44.1
Table 9: Retrieval results measured by nDCG@10 on BEIR (all datasets except MS
MARCO) and the retrieval cluster of MTEB (first 15 datasets). LE stands for
lexical enrichment (our method; §6.2), that enriches query and passage
representation with lexical information.


# Modified Singly-Runge-Kutta-TASE methods for the numerical solution of stiff
differential equations111Partially supported by project PID2022-141385NB-I00
of Ministerio de Ciencia e Innovación of Spain.

M. Calvo<EMAIL_ADDRESS>J. I. Montijano<EMAIL_ADDRESS>L. Rández
<EMAIL_ADDRESS>Departamento Matemática Aplicada, Universidad de Zaragoza.
50009-Zaragoza, Spain.

###### Abstract

Singly-TASE operators for the numerical solution of stiff differential
equations were proposed by Calvo et al. in J.Sci. Comput. 2023 to reduce the
computational cost of Runge-Kutta-TASE (RKTASE) methods when the involved
linear systems are solved by some $LU$ factorization. In this paper we propose
a modification of these methods to improve the efficiency by considering
different TASE operators for each stage of the Runge-Kutta. We prove that the
resulting RKTASE methods are equivalent to $W$-methods (Steihaug and
Wolfbrandt, Mathematics of Computation,1979) and this allows us to obtain the
order conditions of the proposed Modified Singly-RKTASE methods (MSRKTASE)
through the theory developed for the $W$-methods. We construct new MSRKTASE
methods of order two and three and demonstrate their effectiveness through
numerical experiments on both linear and nonlinear stiff systems. The results
show that the MSRKTASE schemes significantly enhance efficiency and accuracy
compared to previous Singly-RKTASE schemes.

###### keywords:

Differential equations; Stiff problems; Runge-Kutta methods; Time-marching
methods; Singly-TASE operators;

##  1  Introduction

Solving stiff initial value problems (IVPs) efficiently and accurately remains
a significant challenge in numerical analysis. Explicit Runge-Kutta (RK)
method, face limitations when applied to stiff systems due to their stability
restrictions. To address these limitations, a new class of time-advancing
schemes for the numerical solution of stiff IVPs was proposed by M. Bassenne,
L. Fu, and A. Mani in [1]. They introduced the concept of TASE (Time-Accurate
and Stable Explicit) operators, designed to enhance the stability of explicit
RK methods. In this approach, instead of solving the differential system

$\dfrac{d}{dt}Y(t)=F(t,Y(t)),\quad Y(t_{0})=Y_{0}\in\mathbb{R}^{d}.$ (1)

they proposed to solve another (stabilized) IVP

$\dfrac{d}{dt}U(t)=T\;F(t,U(t)),\quad U(t_{0})=U_{0}\equiv
Y_{0}\in\mathbb{R}^{d},$ (2)

where $T=T(\Delta t)$ is a linear operator that may depend on the time step
size $\Delta t>0$ and on $F$, so that the numerical solution of Eq. (2)
$U_{RK}(t_{0}+\Delta t)$ obtained with an explicit RK method of order $p$
satisfies $U_{RK}(t_{0}+\Delta t)-Y(t_{0}+\Delta t)=\mathcal{O}(\Delta
t^{p+1}),$ i.e., approximates the local solution of Eq. (1) with order $p$,
and also satisfies some stability requirements that are necessary for solving
stiff systems such as A- or L- stability. In this way, the introduction of the
TASE operator into the original governing equation (1) allows to overcome the
numerical stability restrictions of explicit RK time-advancing methods for
solving stiff systems.

In terms of the accuracy order of the numerical solution of Eq. (2), if the
explicit RK scheme has order $p$ and if the TASE operator satisfies

$T(\Delta t)=I+\mathcal{O}(\Delta t^{q}),$ (3)

it can be easily proved that the numerical solution of the stabilized system
Eq. (2), $U_{RK}(t_{n})$, has at least order $\min(p,q)$. Some specific
schemes combining explicit $p$–stage RK methods with orders $p\leq 4$ and TASE
operators with the same order were derived in [1].

A more general family of TASE operators was derived in [6] by taking

$T(\Delta t)=\sum_{j=1}^{p}\beta_{j}(I-\alpha_{j}\Delta tW)^{-1},$ (4)

where $\alpha_{j}>0$, $j=1,\ldots,p$, are free real parameters, and
$\beta_{j}$ are uniquely determined by the condition $T_{p}(\Delta
t)=I+\mathcal{O}(\Delta t^{p})$. The free parameters $\alpha_{j}$ were
selected to improve the linear stability properties of an explicit RK method
with order $p$ for the scaled equation Eq. (2).

In the TASE operators described by Eq. (4), each evaluation of $T$ involves
the solution of $p$ linear systems with matrices $(I-\Delta
t\;\alpha_{j}\;W)$, $j=1,\ldots,p$. If these systems are to be solved by some
$LU$ matrix factorization, it would be much more efficient if the $\alpha_{j}$
coefficients were all equal, but this is not compatible with order $p$ for the
TASE operator.

To make the methods more efficient, Calvo et al. [2] developed a different
family of TASE operators, termed Singly-TASE operators, defined by

$T_{p}(\Delta t)=\sum_{j=1}^{p}\beta_{j}(I-\alpha\Delta tW)^{-j}.$ (5)

These operators lead to significant improvements in computational efficiency
while maintaining the stability and accuracy necessary for solving stiff
problems.

For an $s$–stage explicit RK method defined by the Butcher tableau

$\begin{tabular}[]{c|ccccc}0&&&&&\\\ $c_{2}$&$a_{21}$&&&&\\\ \vdots&&&&&\\\
$c_{s}$&$a_{s1}$&$a_{s2}$&\ldots&$a_{s,s-1}$&\\\
\hline\cr&$b_{1}$&$b_{2}$&\ldots&$b_{s-1}$&$b_{s}$\end{tabular}\qquad
c_{i}=\sum_{j=1}^{i-1}a_{ij},\quad i=2,\ldots,s,$

the numerical solution of Eq. (2) by a Singly-RKTASE method is advanced from
$(t_{0},U_{0})\to(t_{1}=t_{0}+\Delta t,U_{1})$ by the formula

$U_{1}=U_{0}+\Delta t\left[b_{1}\;K_{1}+\ldots+b_{s}\;K_{s}\right],$

where $K_{j}$, $j=1,\ldots,s$ are computed recursively from the formulas

$\displaystyle K_{1}$ $\displaystyle=$ $\displaystyle T\;F(t_{0},U_{0}),$
$\displaystyle K_{2}$ $\displaystyle=$ $\displaystyle T\;F(t_{0}+c_{2}\Delta
t,U_{0}+\Delta t\;a_{21}K_{1}),$ $\displaystyle\vdots$ $\displaystyle K_{s}$
$\displaystyle=$ $\displaystyle T\;F\left(t_{0}+c_{s}\Delta t,U_{0}+\Delta
t\;\sum_{j=1}^{s-1}a_{sj}K_{j}\right).$ (6)

In this paper, we propose a further modification of the Singly-RKTASE methods
to enhance efficiency, so that for the ith stage of the RK method we use a
TASE operator defined by

$T_{i}=\sum_{j=1}^{r}\beta_{ij}\;(I-\Delta t\;\alpha\;W)^{-j},\quad\alpha>0.$
(7)

The coefficients $\alpha>0$, $r$ and $\beta_{i,j}$, $j=1,\ldots,r$ must be
determined so that the resulting modified singly-RKTASE (MSRKTASE) method

$\begin{array}[]{lll}K_{1}&=&T_{1}\;F(t_{0},U_{0}),\\\
K_{2}&=&T_{2}\;F(t_{0}+c_{2}\Delta t,U_{0}+\Delta t\;a_{21}K_{1}),\\\
\vdots&&\\\ K_{s}&=&T_{s}\;F\left(t_{0}+c_{s}\Delta t,U_{0}+\Delta
t\;\sum_{j=1}^{s-1}a_{sj}K_{j}\right)\\\ U_{1}&=&U_{0}+\Delta
t(b_{1}K_{1}+\ldots+b_{s}K_{s})\end{array}$ (8)

has order $p$. Compared to Singly-RKTASE methods, instead of $s$ coefficients
$\beta_{i}$ we have $sr$ coefficients $\beta_{ij}$, providing more freedom to
improve accuracy and stability properties. Although we could take different
$r_{i}$ for each stage, we have taken the same $r$ for simplicity. Note that
if we impose that all operators $T_{i}$ in (7) have order $p$ with $r=p$, the
coefficients $\beta_{ij}$ are uniquely determined, resulting in a Singly-
RKTASE. Then, we must consider $r>p$, which increases computational cost.
Alternatively, it is possible to use operators $T_{i}$ with order $q<p$ in
such a way that the resulting modified Singly-RKTASE method has order $p$.
This requires further analysis of the conditions the methods’s coefficients
must satisfy to achieve order $p$.

The rest of the paper is organized as follows: In section 2, we prove that
modified Singly-RKTASE methods are a particular case of W-methods. In section
3, we derive the conditions that the coefficients $\beta_{ij}$ and $\alpha$ of
the MSRKTASE must satisfy to have order $p$. In section 4 we develop MSRKTASE
methods of orders 2 and 3 with optimal stability and accuracy properties.
Finally, in section 5, through numerical experiments, we demonstrate the
effectiveness of the newly developed methods, highlighting their potential for
solving both linear and nonlinear stiff systems with improved performance.

##  2 Equivalence between modified Singly-RKTASE and W-methods

Given a matrix $W$ (usually some approximation to the Jacobian matrix at the
current integration point $(t_{n},Y_{n})$) a singly $W$-method [11] with
$\hat{s}$ stages provides an approximation to the solution of Eq. (1) at
$t_{n+1}=t_{n}+h$ by the formulas

$\begin{array}[]{l}(I-h\alpha W)\widehat{K}_{1}=hf(Y_{n})\\\ (I-h\alpha
W)\widehat{K}_{i}=hf(Y_{n}+\sum_{j=1}^{i-1}\hat{a}_{ij}\widehat{K}_{j})+h\alpha
W\sum_{j=1}^{i-1}l_{ij}\widehat{K}_{j},\qquad i=2,\ldots,\hat{s}\\\
Y_{n+1}=Y_{n}+\sum_{i=1}^{\hat{s}}\hat{b}_{i}\widehat{K}_{i}\end{array}$

The coefficients of the method can be arranged in a Butcher tableau as

$\begin{array}[]{c|c|c}\hat{c}&\hat{A}&L\\\
\hline\cr&\hat{b}^{T}&\alpha\end{array}\quad=\quad\begin{array}[]{c|cccc|
ccccc}0&&&&&&\\\ \hat{c}_{2}&\hat{a}_{21}&&&&l_{21}\\\
\vdots&&\ddots&&&\vdots&\ddots\\\
\hat{c}_{s}&\hat{a}_{s1}&\ldots&\hat{a}_{s,s-1}&&l_{s1}&\ldots&l_{s,s-1}\\\
\hline\cr&\hat{b}_{1}&\ldots&\hat{b}_{s-1}&\hat{b}_{s}&\alpha\end{array}$

On the other side, the product of the operator $T_{i}$ in (7) by a vector $v$
can be written as

$\begin{array}[]{l}w_{1}=(I-h\alpha W)^{-1}v\\\ w_{j}=(I-h\alpha
W)w_{j-1},\quad j=2,\ldots,r\\\
T_{i}v=\beta_{i1}w_{1}+\ldots+\beta_{ir}w_{r}\end{array}$

Then, the equations (7)(8) defining a Modified Singly-RKTASE scheme can be
written as

$\begin{array}[]{l}(I-h\alpha W)\widehat{K}_{1}=hf(Y_{n})\\\ (I-h\alpha
W)\widehat{K}_{j}=K_{j-1},\quad j=2,\ldots,r\\\
K_{1}=\beta_{11}\widehat{K}_{1}+\ldots+\beta_{1r}\widehat{K}_{r}\\\ (I-h\alpha
W)\widehat{K}_{r(i-1)+1}=hf(Y_{n}+\sum_{j=1}^{i-1}a_{ij}K_{j})\qquad
i=2,\ldots,s\\\ (I-h\alpha W)\widehat{K}_{r(i-1)+j}=K_{r(i-1)+j-1},\quad
j=2,\ldots,r\\\
K_{i}=\beta_{i1}\widehat{K}_{r(i-1)+1}+\ldots+\beta_{ir}\widehat{K}_{ri}\\\
Y_{n+1}=Y_{n}+\sum_{i=1}^{s}b_{i}K_{i}\end{array}$

Taking into account that $\widehat{K}_{j}=\widehat{K}_{j-1}+h\alpha
W\widehat{K}_{j},$ we have

$\begin{array}[]{l}(I-h\alpha W)\widehat{K}_{1}=hf(Y_{n})\\\ (I-h\alpha
W)\widehat{K}_{j}=hf(Y_{n})+h\alpha W\sum_{l=1}^{j-1}K_{l},\quad
j=2,\ldots,r\\\
K_{1}=\beta_{11}\widehat{K}_{1}+\ldots+\beta_{1r}\widehat{K}_{r}\\\ (I-h\alpha
W)\widehat{K}_{r(i-1)+1}=hf(Y_{n}+\sum_{j=1}^{i-1}a_{ij}K_{j})\qquad
i=2,\ldots,s\\\ (I-h\alpha
W)\widehat{K}_{r(i-1)+j}=hf(Y_{n}+\sum_{j=1}^{i-1}a_{ij}K_{j})+h\alpha
W\sum_{l=1}^{j-1}K_{r(i-1)+l-1},\quad j=2,\ldots,r\\\
K_{i}=\beta_{i1}\widehat{K}_{r(i-1)+1}+\ldots+\beta_{ir}\widehat{K}_{ri}\\\
Y_{n+1}=Y_{n}+\sum_{i=1}^{s}b_{i}K_{i}\end{array}$

and this leads to

$\begin{array}[]{l}(I-h\alpha W)\widehat{K}_{1}=hf(Y_{n})\\\ (I-h\alpha
W)\widehat{K}_{j}=hf(Y_{n})+h\alpha W\sum_{l=1}^{j-1}K_{l},\quad
j=2,\ldots,r\\\ (I-h\alpha
W)\widehat{K}_{r(i-1)+j}=hf(Y_{n}+\sum_{j=1}^{i-1}a_{ij}\sum_{l=1}^{r}\beta_{jl}\widehat{K}_{r(i-1)+j})+h\alpha
W\sum_{l=1}^{j-1}K_{r(i-1)+l-1},\quad i=2,\ldots,s,\quad j=2,\ldots,r\\\
Y_{n+1}=Y_{n}+\sum_{i=1}^{s}b_{i}\sum_{l=1}^{r}\beta_{il}\widehat{K}_{r(i-1)+j}\end{array}$

This is just an $W$ method with $sr$ stages whose vector $\widehat{b}$ and
matrices $\widehat{A}$, $L$ are given by

$L=\begin{pmatrix}L_{r}&&\\\ &\ddots\\\
&&L_{r}\end{pmatrix}\in\mathbb{R}^{rs\times
rs},\qquad\hat{A}=\begin{pmatrix}0&&\\\ A_{21}&0\\\ \vdots&\\\
A_{s1}&\cdots&A_{s,s-1}&0\end{pmatrix}\in\mathbb{R}^{rs\times rs},$
$\widehat{b}=(b_{1}\beta_{11},\ldots,b_{1}\beta_{1r},\ldots,b_{s}\beta_{s1},\ldots,b_{s}\beta_{sr})^{T}$

with

$L_{r}=\begin{pmatrix}0&&\\\ 1&0&\\\ \vdots&\ddots&\ddots\\\
1&\cdots&1&0\end{pmatrix}\in\mathbb{R}^{r\times r},\qquad
A_{i,j}=a_{ij}\begin{pmatrix}\beta_{j1}&\ldots&\beta_{jr}\\\ \vdots&\\\
\beta_{j1}&\ldots&\beta_{jr}\end{pmatrix}\in\mathbb{R}^{r\times r}$

This equivalence let us analyze the order of the Modified Singly-RKTASE
through the order conditions of W-methods. Also, the absolute stability can be
studied through the $W$-methods. The stability function of an $W$ method is
given by

$\hat{R}(z)=1+z\hat{b}^{T}(I-z(\hat{A}+\Gamma))^{-1}\mathbf{1}$

and the limit of this function when $z$ goes to infinity is given by

$\lim_{z\to\infty}\hat{R}(z)=1-\hat{b}^{T}(\hat{A}+\Gamma)^{-1}\mathbf{1}.$
(9)

We will use this to develop Modified Singly-RKTASE methods. Recall that

  * 1.

A method is called A-stable if $|\hat{R}(z)|\leq 1$ for all Re $z\leq 0$;

  * 2.

A method is said to be L–stable if it is A–stable and $\hat{R}(\infty)=0$;

  * 3.

A method is called A$(\theta)$–stable if $|\hat{R}(z)|<1$ for all $z$ such
that $\arg(-z)\leq\theta$, that is, its stability region contains the left
hand region of the complex plane with angle $\theta$. If in addition
$\hat{R}(\infty)=0$, it is called L$(\theta)$–stable.

##  3 Order conditions

Denoting $\Gamma=\alpha(I+L)$, and $\mathbf{1}=(1,\ldots,1)^{T}$, the order
conditions for an W-method are given in Table 1 (see e.g. [7, 11]).

Table 1: order conditions for $W$-methods up to order 4
$\begin{array}[]{|ll|lll|}\hline\cr\hat{b}^{T}\mathbf{1}=1&&&&\\\
\hline\cr\hat{b}^{T}\hat{c}=1/2&&\hat{b}^{T}\Gamma\mathbf{1}=0&&\\\
\hline\cr\hat{b}^{T}\hat{c}^{2}=1/3&\hat{b}^{T}\hat{A}\hat{c}=1/6&\hat{b}^{T}\Gamma^{2}\mathbf{1}=0&\hat{b}^{T}\hat{A}\Gamma\mathbf{1}=0&\hat{b}^{T}\Gamma\hat{A}\mathbf{1}=0\\\
\hline\cr\hat{b}^{T}\hat{c}^{3}=1/4&\hat{b}^{T}\hat{A}\hat{c}\cdot\hat{c}=1/8&\hat{b}^{T}\Gamma^{3}\mathbf{1}=0&\hat{b}^{T}\hat{A}\Gamma^{2}\mathbf{1}=0&\hat{b}^{T}\Gamma\hat{A}\Gamma\mathbf{1}=0\\\
\hat{b}^{T}\hat{A}\hat{c}^{2}=1/12&\hat{b}^{T}\hat{A}^{2}\hat{c}=1/24&b^{T}\Gamma^{2}\hat{A}\mathbf{1}=0&\hat{b}^{T}\hat{A}^{2}\Gamma\mathbf{1}=0&\hat{b}^{T}\hat{A}\Gamma\hat{A}\mathbf{1}=0\\\
&&\hat{b}^{T}\Gamma\hat{A}^{2}\mathbf{1}=0&\hat{b}^{T}\Gamma\hat{c}^{2}=0&\hat{b}^{T}(\hat{A}\Gamma\mathbf{1})\cdot\hat{c}=0\\\
\hline\cr\end{array}$

In the case the operator $W$ is exactly the Jacobian matrix, $W=\partial
f/\partial y(t_{n},Y_{n})$ (Rosenbrock-Wanner methods), many elementary
differentials become the same and the order conditions reduce to those in
Table 2.

Table 2: Order conditions for a Rosenbrock method ($W=\partial f/\partial y$)
up to order 4 $\begin{array}[]{|l|lll|}\hline\cr\hat{b}^{T}\mathbf{1}=1&&&\\\
\hline\cr\hat{b}^{T}(\Gamma+\hat{A})c=1/2&&&\\\
\hline\cr\hat{b}^{T}\hat{c}^{2}=1/3&\hat{b}^{T}(\Gamma+\hat{A})^{2}\mathbf{1}=1/6&&\\\
\hline\cr\hat{b}^{T}\hat{c}^{3}=1/4&\hat{b}^{T}(\hat{A}(\Gamma+\hat{A})\mathbf{1})\cdot\hat{c}=1/8&\hat{b}^{T}(\Gamma+\hat{A})\hat{c}^{2}=1/12&\hat{b}^{T}(\Gamma+\hat{A})^{3}\mathbf{1}=1/24\\\
\hline\cr\end{array}$

Let us analyze the order conditions for a Modified Singly-RKTASE scheme
expressed as its equivalent $W$-method. We can take without losing generality
$\beta_{i1}+\ldots+\beta_{ir}=1$ because this is equivalent to re-scale the
coefficients $b,A$ of the underlying RK method.

First, the order conditions that do not involve the matrix $\Gamma$ reduce to
the order conditions associated to the undelying RK method. For example,
$\hat{b}^{T}\mathbf{1}=b_{1}+\ldots+b_{s}$, or
$\hat{b}^{T}\hat{A}\hat{c}=b^{T}Ac$ which are the corresponding equations of
the RK method.

On the other hand, the form of the matrix $L$ introduces some incompatibility
between some of the order conditions leading to the next

######  Theorem 3.1.

A Modified Singly-RKTASE method can not have order $r+1$.

Proof. Since the matrix $\Gamma$ is a block diagonal matrix with blocks
$\alpha(I_{r}+L_{r})$, it is clear that $(\Gamma-\alpha I)^{r}=(\alpha
L)^{r}=0$ because $L_{r}$ is strictly lower triangular with dimension $r$.
Then, $(\Gamma-\alpha I)^{r}\mathbf{1}=0$. If the method had order $r+1$, by
the order conditions it should be $\widehat{b}^{T}(\Gamma-\alpha
I)^{r}\mathbf{1}=\alpha^{r}$ which is not zero if $\alpha\neq 0$.

With this result, a Modified Singly-RKTASE method with order $p$ must have
necessarily $r\geq p$.

##  4 Modified Singly-RKTASE methods with orders 2 and 3

###  4.1 Methods with $s=2$, $r=2$ and order 2

There exists a family of Runge-Kutta schemes with two stages and order two
that depend on one coefficient $c_{2}$ given by

$a_{21}=c_{2},\quad b_{2}=\dfrac{2-3c3}{6c2(c2-c3)},\quad b_{1}=1-b_{2}$

Since $b^{T}Ac=0$ for any $c_{2}$, the associated order condition of order
three does not depend on $c_{2}$ and the other one vanishes for $c_{2}=2/3$.
So, we will take this value of $c_{2}$.

The only remaining condition for a Modified Singly-RKTASE with $s=r=2$ to
achieve order 2 is $\hat{b}^{T}\Gamma\mathbf{1}=0$ which is satisfied for

$\beta_{22}=-(4+\beta_{12})/3.$

The limit of the stability function when $z$ approaches infinity for the case
of our Modified Singly-RKTASE method of order two is given by

$\lim_{z\to\infty}\hat{R}(z)=1-\hat{b}^{T}(\hat{A}+\Gamma)^{-1}\mathbf{1}=-\dfrac{-7+12\alpha-6\alpha^{2}+6\beta_{12}+\beta_{12}^{2}}{6\alpha^{2}}.$

There are two values of $\beta_{12}$ that make $\hat{R}(\infty)=0$. The one
that yields better results in terms of stability and accuracy is

$\beta_{12}=-3+\sqrt{16-12\alpha+6\alpha^{2}}$

There is only one free parameter $\alpha$ left, which we will use to maximize
the stability region and minimize the error coefficients. The norm of the
error coefficients of order 3 is given by

$C_{3}=\left((\hat{b}^{T}\hat{A}\hat{c}-1/6)^{2}+(\hat{b}^{T}\Gamma^{2}\mathbf{1})^{2}+(\hat{b}^{T}\hat{A}\Gamma\mathbf{1})^{2}+(\hat{b}^{T}\Gamma\hat{A}\mathbf{1})^{2}\right)^{1/2}.$

The corresponding norm when $W=f^{\prime}$ is

$D_{3}=\left|\hat{b}^{T}(\Gamma+\hat{A})^{2}\mathbf{1}-1/6\right|.$

The method is $A$-stable (and therefore $L$-stable) if
$\alpha\in[0.3117,3.257]$ and $C_{3}$ is monotonic decreasing with $\alpha$. A
good value of the parameter can is $\alpha=32/100$. With these values, we
obtain an $L$-stable method of order 2, which we will name MSRKTASE2. The
error coefficients are given in Table 3. For comparison, we also include the
error coefficients of the Singly-RKTASE of order 2 obtained in [2], named
SRKTASE2. As we can see, the new method achieves $L$-stability (SRKTASE2 has
$\hat{R}(\infty)=1/2$) and error coefficients that are about 20 times smaller
(40 times smaller in the case where $W=f^{\prime}$).

We could use the parameter $\beta_{12}$ to obtain smaller error coefficients
without imposing the condition $R(\infty)=0$. For example we could nullify the
coefficient $D_{3}$ (order three for the case $W=f^{\prime}$,) but in this
case, we only achieve reasonable stability regions for large values of
$\alpha$ and $C_{3}$ becomes very large (greater than one hundred).

Table 3: properties of the proposed methods Method | $p$ | $\alpha$ | $C_{p+1}$ |  $D_{p+1}$ ( $W=f^{\prime}$ ) |  $\theta$
---|---|---|---|---|---
SRKTASE2 | $2$ | $2$ | $4.00347$ | $4.16667$ |  $90^{\circ}$
MSRKTASE2 | $2$ | $0.32$ | $0.212866$ | $0.10116$ |  $90^{\circ}$
SRKTASE3 | $3$ | $1.8868$ | $6.7171$ | $6.6753$ |  $88.99^{\circ}$
MSRKTASE3a | $3$ | $0.54$ | $0.1817$ | $0.2288$ |  $88.23^{\circ}$
MSRKTASE3b | $3$ | $0.56$ | $0.3968$ | $0.0035$ | $50.38^{\circ}$

###  4.2 Methods with $s=3$, $r=3$ and order 3

There exist a family of Runge-Kutta schemes with three stages and order three,
defined by the coefficients

$a_{32}=\dfrac{c_{3}(c_{2}-c_{3})}{c_{2}(3c_{2}-2)},\quad
b_{2}=\dfrac{2-3c3}{6c2(c2-c3)},\quad b3=\dfrac{2-3c2}{6c3(c3-c2)},\quad
b_{1}=1-b_{2}-b_{3}$

depending on two free parameters $c_{2}$ and $c_{3}$ with $c_{2}\neq c_{3}$
and $c_{2}\neq 2/3$.

The remaining conditions for a Modified Singly-RKTASE with $s=r=3$ to achieve
order 3 form a set of four equations $\hat{b}^{T}\Gamma\mathbf{1}=0$,
$\hat{b}^{T}\Gamma^{2}\mathbf{1}=0$, $\hat{b}^{T}\hat{A}\Gamma\mathbf{1}=0$,
$\hat{b}^{T}\Gamma\hat{A}\mathbf{1}=0$ which can be solved for the parameters
$\beta_{12},\beta_{13},\beta_{23}$ and $\beta_{33}$, resulting in:

$\begin{array}[]{l}\beta_{12}=\dfrac{c_{3}(-2+3c_{3})\beta_{22}-3c_{2}^{2}(6c_{3}+\beta_{32})+2c_{2}(9c_{3}^{2}+\beta_{32})}{(c_{2}-c_{3})(2-3c_{3}+c_{2}(-3+6c_{3}))},\\\\[9.0pt]
\beta_{13}=-\dfrac{c_{3}(-2+3c_{3})(1+\beta_{22})-3c_{2}^{2}(1+4c_{3}+\beta_{32})+2c_{2}(1+6c_{3}^{2}+\beta_{32})}{2(c_{2}-c_{3})(2-3c_{3}+c_{2}(-3+6c_{3}))}\\\\[7.0pt]
\beta_{23}=1/2(-1-\beta_{22}),\\\ \beta_{33}=1/2(-1-\beta_{32}),\end{array}$

depending on the parameters $c_{2},c_{3},\beta_{22}$ and $\beta_{32}$. Note
that in this family of methods of order three, $\alpha$ is also a free
parameter. Imposing $\beta_{32}=\beta_{22}=\beta_{12}$ we obtain the third
order Singli-RKTASE methods obtained in [2].

The free parameters can be selected to maximize the stability region and to
minimize the coefficients of the leading term of the local error. For these
methods of order three, it is also satisfied that

$\hat{b}^{T}\hat{A}\Gamma\hat{A}\mathbf{1}=0,\quad\hat{b}^{T}\hat{A}^{2}\Gamma\mathbf{1}=0,\quad\hat{b}^{T}\Gamma\hat{A}^{2}\mathbf{1}=0,\quad\hat{b}^{T}\Gamma\hat{A}\Gamma\mathbf{1}=0,\quad\hat{b}^{T}\Gamma\hat{c}^{2}=0,\quad\hat{b}^{T}(\hat{A}\Gamma\mathbf{1})\cdot\hat{c}=0$

and the other error coefficients of the term of order four are

$\begin{array}[]{l}C_{41}\equiv\hat{b}^{T}\hat{A}^{2}\hat{c}-1/24=b^{T}A^{2}c-1/24=-1/24\\\\[5.0pt]
C_{42}\equiv(8\hat{b}^{T}\hat{A}\hat{c}\cdot\hat{c}-1)/24=(8b^{T}Ac\cdot
c-1)/24=(4c_{3}-3)/72\\\\[5.0pt]
C_{43}\equiv(12\hat{b}^{T}\hat{A}\hat{c}^{2}-1)/24=(12b^{T}Ac^{2}-1)/24=(2c_{2}-1)/24\\\\[5.0pt]
C_{44}\equiv(4\hat{b}^{T}\hat{c}^{3}-1)/24=(4b^{T}c^{3}-1)/24=(2((c_{2}(2-3c_{3})+2c_{3})-3)/72\\\\[6.0pt]
C_{45}\equiv\hat{b}^{T}\hat{A}\Gamma^{2}\mathbf{1}=\alpha^{2}\frac{-2\beta_{22}+6c_{3}(3+\beta_{22})-3c_{3}^{2}(3+\beta_{22})+2\beta_{32}+9c_{2}^{2}(3+\beta_{32})-3c_{2}(6-\beta_{22}+2c_{3}(3+\beta_{22})+3\beta_{32}))}{12(c_{2}-c_{3})(2-3c_{3}+c_{2}(-3+6c_{3}))}\\\\[6.0pt]
C_{46}\equiv\hat{b}^{T}\Gamma^{2}\hat{c}=\alpha^{2}\frac{(-2\beta_{22}+3c_{3}(3+\beta_{22})+2\beta_{32}-3c_{2}(3+\beta_{32}))}{12(c_{2}-c_{3})}\\\\[7.0pt]
C_{47}\equiv\hat{b}^{T}\Gamma^{3}\mathbf{1}=\alpha^{3}\end{array}$ (10)

If we solve $\hat{b}^{T}\Gamma^{2}\hat{c}=0$ and
$\hat{b}^{T}\hat{A}\Gamma^{2}\mathbf{1}=0$ for $\beta_{22}$ and $\beta_{32}$
we obtain the method in [2], independent of the values of $c_{2}$ and $c_{3}$.

The first four equations in (10) depend only on $c_{2}$ and $c_{3}$ and can
not be satisfied simultaneously. The 2-norm of these four equations reaches
its minimum value at $c_{2}=0.496188$, $c_{3}=0.764887$, very close to
$c_{2}=1/2$, $c_{3}=3/4$ for which

$\dfrac{\big{(}(24b^{T}A^{2}c-1)^{2}+(12b^{T}Ac^{2}-1)^{2}+(8b^{T}Ac\cdot
c-1)^{2}+(4b^{T}c^{3}-1)^{2}\big{)}^{1/2}}{24}=\dfrac{\sqrt{145}}{288}=0.041811$

From now on, we will take $c_{2}=1/2$ and $c_{3}=3/4$. Thus, we have three
coefficients $\alpha,\beta_{22},\beta_{32}$ to minimize the error coefficients
and maximize the stability region.

The limit of the stability function of this method as $z$ approaches infinity
is given by

$\lim_{z\to\infty}\hat{R}(z)=1-\hat{b}^{T}(\hat{A}+\Gamma)^{-1}\mathbf{1}=\dfrac{a_{0}+a_{1}\alpha-288\alpha^{2}+96\alpha^{3}}{96\alpha^{3}}$

with

$\begin{array}[]{l}a_{0}=-(-3+\beta_{22})(-3+\beta_{32})(33+3\beta_{22}+4\beta_{32})\\\
a_{1}=-6(-45+12\beta_{22}+\beta_{22}^{2})\end{array}$

There are two values of $\beta_{32}$ (as functions of $\alpha$ and
$\beta_{22}$) that make $R(\infty)=0$. We will select one of these values and
use the coefficients $\alpha$ and $\beta_{22}$ to minimize the error
coefficients and maximize the $L(\theta)$-stability angle.

The 2-norm of the error coefficients is given by

$C_{4}=(C_{41}^{2}+\cdots+C_{47}^{2})^{1/2}$

Minimizing this norm is equivalent to minimizing (the other coefficients are
constant) $C_{45}^{2}+C_{46}^{2}+C_{47}^{2}.$ A good compromise between large
$\theta$ and small $C_{4}$ is obtained for

$\alpha=0.54,\quad\beta_{22}=-6.1,\quad\beta_{32}=-2.75034$

resulting in a method, that we will name MSRKTASE3a, with stability angle and
error coefficients given in Table 3.

Alternatively, we can minimize the 2-norm of the coefficients of the leading
term in the case that $W$ is exactly the Jacobian matrix,

$D_{4}=(D_{41}^{2}+\cdots+D_{45}^{2})^{1/2}.$

This is equivalent to minimize
$|D_{41}|=\hat{b}^{T}(\Gamma+\hat{A})^{3}\mathbf{1}-1/24$ (the other
coefficients are constant). This error coefficient can be made to vanish when

$\beta_{22}=-3-\dfrac{1}{3\alpha^{2}}+8\alpha.$

A good compromise between maximizing $\theta$ and minimizing $C_{4}$ with
$\alpha$ is obtained when

$\alpha=0.56,\quad\beta_{22}=-6.1,\quad\beta_{32}=-2.75034$

resulting in a method, that we will name MSRKTASE3b, with stability angle and
error coefficients given in Table 3.

In Figure 1 we plot the boundary of the stability regions of the proposed
methods. The method with minimal $D_{4}$ is shown in red-dashed, the method
with minimal $C_{4}$ in red-solid, and the method from [2] in blue.

Figure 1: Stability regions of the new methods (left). On the right, a zoom of
the region near the origin.

We can see from these results that the new method MSRKTASE3a, has similar
stability properties to the method of order 3 in [2], but has about 35 times
smaller error coefficients. The new method is expected to provide more
accurate approximations with the new method. The other method, MSRKTASE3b,
also has much smaller error coefficients, although the stability angle is
smaller. Nonetheless, its stability region is not significantly worse.

##  5 Numerical experiments

To evaluate the performance of the new methods, we considered two test
problems and integrated them using the two new methods, MSRKTASE3a and
MSTASE3b. For comparison, we also used the method proposed in [2] (STASE) to
demonstrate the improvements of these new modified singly-RKTASE methods over
the previous singly-RKTASE methods. Additionally, to benchmark the performance
of the new methods against other known Runge-Kutta methods for stiff problems,
we integrated the problems with a Singly-Diagonally Implicit RK method of
order 3 (SDIRK) proposed in [8]. Since the SDIRK method is implicit, we used a
simplified Newton method to solve the stage equations, approximating the
Jacobian matrix with the matrix $W$ used in the TASE operators.

For each method and problem, we computed the CPU time required for solving it
and the 2-norm of the global error at the end of the integration interval.
These data points allowed us to plot the global error (in logarithmic scale)
against the CPU time, resulting in an efficiency plot.

Problem 1: (taken from [1] and [6]) The 1D diffusion of a scalar function
$y=y(x,t)$, with a time dependent source term

$\dfrac{\partial y}{\partial t}=\dfrac{\partial^{2}y}{\partial
x^{2}}+0.1\sin(t/50),\qquad y(x,0)=1-\cos(x)^{101},\qquad 0\leq x\leq 2\pi.$

The solution $y=y(x,t)$ is assumed to be $2\pi$-periodic in $x$. For the
spatial discretization, we used fourth-order centered difference schemes with
a grid resolution of $N=512$, assuming periodic boundary conditions. The real
part of the eigenvalues of the Jacobian matrix of the semi-discrete problem
ranges from $-3.5\times 10^{4}$ to $2.79\times 10^{-5}$.

The efficiency plot of the methods for this problem is depicted in Figure 2.

Figure 2: Efficiency plot for Diffusion problem

The method MSRKTASE3b is as efficient as SDIRK and more efficient than
MSRKTASE3a. This is due to the fact that the eigenvalues of the Jacobian
matrix are close to the real axis, therefore included in the stability region
and the Jacobian matrix $J$ of the problem is constant, then $W=J$. In this
situation MSRKTASE3b has smaller error coefficients than MSRKTASE3a and this
translates into a better efficiency. MSRKTASE3a is clearly more efficient than
SRKTASE.

Problem 2: The 1D Burgers’ equation written in the conservative form

$\dfrac{\partial y}{\partial t}=0.1\dfrac{\partial^{2}y}{\partial
x^{2}}-\dfrac{\partial}{\partial x}\left(\dfrac{y}{2}\right)^{2},\qquad
y(x,0)=1-\cos(x)^{101},\qquad 0\leq x\leq 2\pi.$

The solution $y=y(x,t)$ is assumed to be $2\pi$-periodic in $x$. For the
spatial discretization, a fourth-order centered difference scheme with grid
resolution of $N=512$.

The real part of the eigenvalues of the Jacobian matrix of the semi-discrete
problem ranges from $-3.5\times 10^{3}$ to $2.6\times 10^{-6}$.

We integrated Burger’s problem using the three methods with three different
options for the matrix $W$. Firstly, the Jacobian matrix is evaluated at every
time step and $W=\partial f(t_{n},y_{n})/\partial y$. With this option, the
$LU$ matrix factorization had to be computed at every step, increasing the
computational cost. Secondly, the Jacobian matrix was evaluated only at the
initial time step and $W=\partial f(t_{0},y_{0})/\partial y$, requiring only
one $LU$ factorization and considerably reducing the CPU time. However, this
could affect the accuracy in the case of TASE methods or the number of
iterations required to solve the non linear system in the case of the SDIRK
method. Finally, $W$ was taken as the linear part of the semi-discrete
differential equation (the matrix of the diffusion term). The computational
cost is similar to the second option, but the accuracy could be lower. This
option provided insight into the methods behaviour when $W$ poorly
approximates the Jacobian matrix.

Figure 3: Efficiency plot for Burger’s problem integrated with $W$ the
Jacobian matrix evaluated at every step (left) and $W$ the Jacobian matrix
evaluated just at the initial time step (right)

Figures 3 and 4 show the efficiency plots for Burger’s problem integrated with
these three options.

As shown in the Figures, evaluating the Jacobian matrix only at the initial
step reduces the CPU time but increases the global error for the TASE-based
methods, where the local error depends on $W$ and is smaller when it coincides
with the Jacobian matrix. For the SDIRK method, the error does no vary
significantly with different $W$, but solving the nonlinear system requires
more iterations when the Jacobian matrix is less accurately approximated,
thereby increasing the computational cost.

Figure 4: Efficiency plot for Burger’s problem integrated with $W$ equal to
the matrix of diffusion term

MSRKTASE3a exhibited stability issues with the two largest time steps when $W$
is not the Jacobian matrix at every step, likely due to a high sensitivity of
the stability function to parameter changes. MSRKTASE3b also showed stability
problems with the largest time step size and the poorest approximation of the
Jacobian matrix. Further research is being carried out in this regard.

Concerning the efficiency of the methods, when the jacobian matrix is
evaluated at every step, SDIRK3 is the most efficient, followed by MSRKTASE3b
and MSRKTASE3a. SRKTASE3 is the least efficient due to its larger error
coefficients. When the Jacobian matrix is evaluated only at the initial point,
MSRKTASE3a and MSRKTASE3b are more efficient while SRKTASE3 remains the least
efficient. Finally, when the matrix $W$ is the matrix of the diffusion term,
MSRKTASE3a is more efficient for small time steps, where it has no stability
problems, due to the relevance of the error coefficient C4. MSRKTASE3b remains
more efficient than SDIRK and SRKTASE3 is also more efficient than SDIRK
except for the smallest time step. The solution of the nonlinear equations
makes SDIRK less efficient.

##  6 Conclusions

We presented a modification of the Singly-RKTASE methods for the numerical
solution of stiff differential equations by taking different TASE operators
for each stage of the RK method. We proved that these methods are equivalent
to $W$-methods which enabled us to derive the order conditions when the TASE
operators have order smaller than the order $p$ of the Runge-Kutta scheme. For
the case where $p=3$ we obtained Modified Singly-RKTASE methods that
significantly reduce the error coefficients compared to Singly-RKTASE methods,
while maintaining stability.

Numerical experiments on both linear and nonlinear stiff systems demonstrate
that the modified Singly-RKTASE methods provide significant improvements in
accuracy and computational efficiency over the original Singly-RKTASE schemes,
making them competitive with other methods such as Diagonally Implicit RK
methods.

## References

  * [1] M. Bassenne, L. Fu, and A. Mani. Time-Accurate and highly-Stable Explicit operators for stiff differential equations. Journal of Computational Physics 424 (2021): 109847.
  * [2] M. Calvo, L. Fu, J. I. Montijano and L. Rández. Singly TASE operators for the numerical solution of stiff differential equations by explicit Runge-Kutta schemes. Journal of Scientific computing (2023)
  * [3] K. Burrage. A special family of RK methods for solving stiff differential equations. BIT 18, 1 (1978): 22-24.
  * [4] J.C. Butcher. On the implementation of implicit RK methods. BIT 16, 3 (1976): 237–240.
  * [5] J. C. Butcher. Numerical methods for ordinary differential equations. John Wiley & Sons, 2003.
  * [6] M. Calvo, J.I. Montijano, L. Rández. A note on the stability of time-accurate and highly-stable explicit operators for stiff differential equations. Journal of Computational Physics 436 (2021): 110316.
  * [7] G. Wanner, E. Hairer. Solving ordinary differential equations II. Springer Berlin Heidelberg, 1996.
  * [8] C.A. Kennedy, M.H. Carpenter. Diagonally Implicit Runge–Kutta Methods for Ordinary Differential Equations: A Review. NASA/TM-2016-219173 NASA Langley Research Center (2016): 1–162
  * [9] C.A. Kennedy, M.H. Carpenter. Diagonally implicit Runge–Kutta methods for stiff ODEs, Applied Numerical Mathematics 146 (2019): 221–244.
  * [10] H.H. Rosenbrock. Some general implicit processes for the numerical solution of differential equations. The Computer Journal 5 (1963): 329–330.
  * [11] T. Steihaug and A. Wolfbrandt. An attempt to avoid exact Jacobian and nonlinear equations in the numerical solution of stiff differential equations. Mathematics of Computation 33 (1979): 521–534.